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ABSTRACT 

This  paper  presents  a  method  for  finding  the  state  of  the  detonation 
products  for  a  plane  detonation  wave  utilizing  the  two-dimensional  steady- 
state  hydrodynamics  equations  and  the  method  of  characteristics.  It  is  shown 
that  in  those  cases  where  the  equation  of  state  gives  pressure  along  the  adiabat 
as  a  function  of  density  (only),  the  characteristic  equations  can  be  solved 
analytically  in  velocity  space  for  the  plane  problem.  A  numerical  method  for 
the  solution  in  coordinate  space  is  then  developed  using  the  solution  in  velocity 
space  in  conjunction  with  the  geometric  relationship  between  characteristic 
curves  in  the  two  spaces.  In  this  manner  the  difficulty  encountered  due  to  the 
coincidence  of  the  detonation  front  and  two  characteristic  curves  in  coordinate 
space  is  overcome.  The  derivations  involved  are  included  as  appendices. 

The  body  of  the  paper  describes  the  application  of  the  technique  for  an  ideal 
gas  equation  of  state  and  for  an  equation  of  state  developed  by  Mark  Wilkins 
(UCRL-7797).  A  further  application  is  described  in  which  this  solution  is 
used  to  start  a  solution  for  a  steady  state  detonation  in  cylindrical  geometry. 
The  results  are  reported  for  PBX  9404  in  the  form  of  graphs  produced  by  the 
computing  facilities  on  a  cathode  ray  tube  from  the  computer  programs 
described. 


GENERAL  INTRODUCTION 

In  the  course  of  work  at  the  Lawrence  Radiation  Laboratory  (Livermore) 
with  the  numerical  calculation  of  the  detonation  of  high  explosives  employing 
finite  difference  techniques  directly  on  the  differential  equations  of  continuum 
mechanics  as  described  in  reference  1  (the  HEMP  code),  it  was  decided  that 
some  kind  of  independent  check  on  the  ability  of  HEMP  to  calculate  a  steady 
state  detonation  was  desirable.  The  objective  was  to  establish  that  HEMP 
could  sustain  a  steady  state  solution  over  a  "long"  time  and  to  establish 
confidence  in  the  numerical  values  of  the  variables  calculated.  It  is  well 
known  that  the  method  of  characteristics  enjoys  the  confidence  of  a  large 
number  of  people  and  would,  therefore,  be  a  good  method  to  accomplish  both 
objectives.  The  work  described  here  was  undertaken  with  this  in  mind. 

It  will  be  noted  that  the  computer  programs  written  lack  "polish"  and 
cannot  be  used  in  present  form  as  production  codes  to  generate  solutions  to 
all  the  "practical"  problems  one  would  desire.  This  lack  of  polish  should  be 
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judged  in  the  light  of  the  objectives  described  in  the  paragraph  above.  The 
work  was  directed  toward  demonstrating  that  the  method  could  be  employed, 
and  then  working  a  representative  problem  which  could  be  used  to  check 
against  the  same  problem  solved  on  HEMP.  These  objectives  have  been  met. 
(This  comparison  will  be  made  in  a  report  by  Mark  Wilkins  to  be  issued  later.) 

The  report  itself  has  been  written  with  an  eye  toward  self- containment. 

It  seemed  desirable  to  include  everything  that  is  necessary  for  a  thorough 
understanding  of  the  development  of  the  numerical  methods  employed.  For 
the  most  part  this  work  has  been  confined  to  the  appendices. 

Sections  I  and  II  are  devoted  to  the  solution  of  a  steady  state  detonation 
in  plane  geometry  and  section  III  describes  the  steady  state  problem  in 
cylindrical  geometry.  The  check  of  HEMP  was  to  be  made  with  the  plane 
problem,  and  the  cylindrical  problem  was  included  because  it  appeared  that 
a  need  was  developing  for  a  production  code  for  cylindrical  steady- state 
detonations  in  connection  with  experimental  programs  under  way  at  the 
Lawrence  Radiation  Laboratory  studying  material  properties.  Such  a  code 
employing  the  method  of  characteristics  would  be  considerably  faster  than 
HEMP  for  this  problem.  The  program  described  merely  demonstrates  that 
the  technique  of  section  III  could  be  employed  to  build  such  a  code. 

I  take  this  opportunity  to  offer  my  thanks  to  Mark  Wilkins  for  proposing 
the  project  and  encouraging  me  throughout,  to  John  Hardy  for  his  suggestions 
that  helped  transform  my  thoughts  into  deeds,  to  Fred  Fritsch,  Don  Emery, 
and  especially  Gloria  Scoggin  for  their  instructions  and  aid  in  programming 
and  debugging,  and  to  the  many  people  of  the  Technical  Information  Division 
for  their  talent  and  patience. 

GENERAL  DISCUSSION 

It  is  generally  accepted  that  the  state  of  the  detonation  products  of  a 
high  explosive  (H.  E.  )  can  be  accurately  calculated  using  the  hydrodynamics 
equations  behind  the  detonation  front  and  the  Hugoniot  equations  across  the 
front  in  conjunction  with  the  Chapman- Jouguet  hypothesis  (see  ref.  2).  It  is 
well  known  that  the  flow  behind  the  front,  relative  to  the  front,  is  supersonic, 
and  that  if  the  detonation  products  are  assumed  to  be  nonviscous  and  non¬ 
conducting  while  the  flow  is  assumed  to  be  irrotational  and  isentropic,  then 
the  steady- state  hydrodynamics  equations  are  hyperbolic.  For  hyperbolic 
equations,  the  method  of  characteristics  can  be  employed  in  those  cases 
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where  the  boundary  conditions  are  given  on  a  curve  that  is  not  itself  a 
characteristic  curve.  For  the  two-dimensional  steady- state  detonation  it 
happens  that  the  detonation  front  is  both  a  boundary  curve  and  a  characteristic 
curve.  In  fact  the  front  occurs  where  two  characteristic  curves  become  co¬ 
incident.  Therefore,  to  use  the  method  of  characteristics  for  this  problem, 
it  is  necessary  to  devise  a  procedure  to  generate  information  behind  this 
front  where  the  method  can  be  used. 

Pack  and  Hill  devised  a  method  employing  an  expansion  in  power  series 
from  the  front  (see  ref.  3).  The  method  described  on  the  following  pages 
employs  the  hodograph  transformation  which  transforms  the  detonation  front 
from  a  line  in  coordinate  space  to  a  point  and  a  characteristic  curve  in 
velocity  space  (see  appendix  6  for  details).  It  is  then  possible  to  use  the 
geometric  relation  between  characteristic  curves  in  the  two  spaces  to  solve 
the  problem  numerically.  This  method  requires  an  equation  of  state  for  the 
detonation  products  which  results,  through  Bernoulli' s  equation,  in  a  sound 
speed,  or  relative  volume,  that  is  a  function  of  the  two  velocity  components. 
This  is  the  same  as  requiring  that  the  pressure  be  a  function  of  density,  only, 
along  an  adiabat.  It  was  not  possible  to  make  a  direct  comparison  with  the 
Pack  and  Hill  results  because  the  equation  of  state  they  used  was  tabular  but 
the  "shapes"  of  the  curves  they  give  are  the  same  as  those  shown  in  this 
report.  The  method  used  here  seems  to  be  easier  to  use  and  does  not  in¬ 
volve  as  many  assumptions. 

To  work  the  plane  problem  using  Wilkins'  equation  of  state  (see  ref.  4) 
it  was  necessary  to  evaluate  an  improper  integral  (see  part  B,  section  II) 
that  is  shown  to  converge  in  appendix  10.  As  a  further  check  on  this  integral 
it  was  evaluated  for  an  ideal  gas  using  the  input  constants  shown  in  part  D, 
section  III,  and  the  problem  discussed  in  section  I  was  worked  using  the 
program  described  in  section  II.  The  results  of  this  are  not  included  but 
when  the  plots  obtained  were  compared  with  the  plots  in  section  I,  they  were 
identical.  The  ideal  gas  equation  of  state  can,  therefore,  be  included  in  a 
program  written  for  the  Wilkins  equation  of  state,  and  this  was  done  for  the 
cylindrical  program. 

The  data  displayed  in  Figs.  2-15  and  17-19  were  included  to  illustrate 
the  kind  of  information  available  directly  from  the  programs  as  written. 

With  the  exception  of  Fig.  19,  the  data  were  plotted  by  the  computing  facilities 
at  Lawrence  Radiation  Laboratory  (Livermore)  from  the  tapes  generated  on 
the  CDC  36  00  high  speed  computer  directly  from  the  programs  described. 
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CONCLUSIONS 

The  work  herein  described  demonstrates  that  the  method  of  character¬ 
istics  can  be  employed  to  determine  the  state  of  the  products  of  a  plane 
steady- state  detonation  by  using  the  hodograph  transformation  and  the 
geometric  relationship  between  the  characteristic  curves  in  the  velocity 
space  and  those  in  the  coordinate  space  to  generate  information  behind  the 
detonation  front.  It  is  also  demonstrated  that  with  some  experimentation  it 
would  be  possible  to  use  the  solution  to  the  plane  problem  to  build  a  production 
code  for  "rapid"  solutions  to  cylindrical  steady-state  detonation  problems. 
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I.  A  PLANE  TWO-DIMENSIONAL  DETONATION  FOR  AN  EXPLOSIVE 
WITH  AN  IDEAL  GAS  EQUATION  OF  STATE 

Introduction 

In  parts  A  through  D  below  is  presented  the  solution  of  a  plane  two- 
dimensional  steady- state  detonation  of  an  explosive  with  an  ideal  gas  equation 
of  state  employing  the  method  of  characteristics.  Part  A  gives  an  outline  of 
the  derivation  of  the  basic  equations  governing  the  motion;  part  B  gives  a 
description  of  the  numerical  method  used;  part  C  describes  the  example 
problem  solved;  and  part  D  gives  a  description  of  the  program  card  decks 
available. 


A.  Calculation  of  the  Two-Dimensional  Steady-State  Detonation 
for  an  Ideal-Gas-Type  Explosive 

The  general  equations  governing  the  motion  of  a  steady  state,  non- 
viscous,  nonconducting  medium  in  the  absence  of  body  forces  in  a  plane  are 
(see  appendix  1): 

Conservation  of  Mass 


9(p  u) 

Bx 


o 

ay 


Conservation  of  Momentum 


3p_  9u  ,  9u 
"  ^  "  pu  ^  +  pv 


9p  9v  ,  9v 
~  pU  Fx  +  pV  9y 


Conservation  of  Energy 


o  /  2,2 

3  (  u  +  v 
Pu  -537  (e  + 


■)  +pv^(e 


2  2 
u  +  v 


9  (pu)  3(pv) 
8x  By 


Entropy  Principle 
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where 

p  =  mass  density, 
u  =  velocity  in  the  x  direction, 
v  =  velocity  in  the  y  direction, 

(x,  y)  =  rectangular  coordinates, 
p  =  hydrostatic  pressure, 
e  =  internal  energy  per  unit  mass, 
n  =  entropy  per  unit  mass. 

If  it  is  assumed  that  the  material  has  the  ideal  gas  equation  of  state, 


so  that  only  the  variables  x,  y,  u,  v  are  involved  in  these  equations.  From 
the-  ideal  gas  equation  of  state  it  can  be  shown  that: 
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e  =  -- r  ,  (see  Ref.  2) 

7  -  1 

E  =  Po€ 


(3) 


where 

e  =  internal  energy  per  unit  mass, 
p  pi’  ai  are  values  at  some  point, 

V  =  relative  volume, 
pQ  =  reference  density. 

The  problem  is  then  reduced  to  finding  u  and  v  and  using  the  relations 
above  to  find  the  other  variables  of  interest.  In  this  treatment  the  equations 
(1)  were  solved  numerically  by  the  use  of  the  method  of  characteristics.  The 
equations  can  be  shown  to  be  hyperbolic  (see  appendix  2)  and  if  the  so-called 
hodograph  transformation  is  made,  i.  e.,  if  the  equations  are  recast  with  u,  v 
as  the  independent  variables,  four  equations  result  (see  appendix  3).  Two 
equations  are  for  characteristics  in  the  u,  v,  plane  (r)  and  two  are  for 
characteristics  in  the  x,  y  plane  (C),  which  last  are  also  "compatibility  rela¬ 
tions"  that  hold  along  r  curves.  These  four  relations  are: 


uv  +  a 


along  which 


/  2  2  2 

uv  +  aVu  +  v  -  a  dx  _  dy 

2  2  do-  do- 

u  -  a 


r 


2: 


dv 

du 


avu  +  v 


v 


2 


2 

a 


(4) 


uv 


a Vu2  +  v^  -  a^  dx  _  dy 

~2  2  dp  "  dj3 

u  “a 


along  which 
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where  a  and  /3  are  parameters  along  r.  and  F _  respectively,  provided  that 

2  ,  2  2  .  n 
u  +  v  -  a  >0. 

Now  if  one  lets 

ip+  =  angle  the  curve  makes  with  the  positive  u  axis, 
i p~  =  angle  the  r 2  curve  makes  with  the  positive  u  axis, 

then  one  can  show  (see  appendix  5)  that  the  equations  (4)  may  be  written  as 

r+:  v  =  (tan  ib+)  u  , 
a  a 


“  tan  4j 
Vp  =  (tan  \jj~)  Up  , 


c  •  +  xr  1 

p  tan  ip  p 

pi  A  -J-  _ 

where  A  =  -rr--  ,  and  the  r  and  C  curves  can  be  designated  r,  =  r  ,  C0  =  C  , 

Off  14 

etc.  because  of  the  particular  way  in  which  they  are  related  (see  below). 

-f-  _  _  -}“ 

Note  that  the  C  and  F  curves  are  perpendicular  as  are  C  and  r  . 

The  coordinate  system  used  is  one  in  which  the  detonation  front  is 
considered  to  be  at  rest  and  the  initial  values  of  the  dependent  variables  are 
those  determined  by  the  Chapman- Jouguet  hypothesis  (see  appendix  4).  Only 
the  "upper"  half  of  the  slab  is  considered  since  the  other  half  is  the  reflection 
of  the  one  considered. 

In  appendices  6  and  7  it  is  shown  that  the  equations  for  r  and  r  can  be 

solved  in  parametric  form.  They  are  shown  to  be  epicycloids  generated  by  a 

circle  of  radius  q  (see  below)  rolling  on  a  circle  of  radius  a  .,  where 

c  ] 

a  2  the  Chapman- Jouguet  sound  speed, 

2  7  -  1 

y  =  specific  heat  ratio. 

The  T+  curves  are  generated  if  the  circle  of  radius  q  rotates  counterclock¬ 
wise,  and  the  F-  curves  are  generated  if  it  rolls  clockwise.  The  constant  in 
equation  (2)  is  determined  in  appendix  6  and  designated  q  =  a  ./g,  and  is 

C-] 
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shown  to  be  the  locus  of  points  where  the  sound  speed  is  zero.  There  are  two 

parameters  involved  in  the  solutions.  The  first  is  the  angle  subtended  at  the 

center  of  the  circle  of  radius  a  .  by  the  center  of  the  rolling  circle  when  it 

^  3 

is  in  its  initial  position  (ip  J.  The  second  is  the  angle  subtended  at  the  center 

of  the  circle  of  radius  a  .  by  the  center  of  the  rolling  circle  as  it  rolls  (/3), 

c  3 

measured  with  respect  to  the  initial  position.  In  terms  of  these  parameters 
the  r  curves  are  given  by 

u  =  (a  .  +  q)  cos(/3  +  i //>:<)  -  q 
C~3 


COS 


a  .  +  q  > 

-^—0  +  0* 


v  =  (a  .  +  q)  sin(/3  +  (//.,.)  -  q  sin 
c-j 


'a  •  +  q 
c-3 

q 


(6) 


0  + 


For  (//.,,  <0,  ft  >  0  the  r+  curves  are  generated,  and  for  >  0,  j3  <  0  the  r 
curves  are  generated. 

In  terms  of  these  same  variables  the  C  curves  are  given  by 


'0 


tan 


0 


1  -  q 

.+ 


X 


0  * 


(7) 


For  <//0,  >  0,  /3  <  0  the  C  curves  are  generated,  and  for  i //.,.<  0,  j3  >  0  the  C 
curves  are  generated. 

It  can  be  shown  (see  appendix  6)  that  for  this  type  of  problem  the  solu- 

2  2  A2 

tion  in  the  u,  v  plane  is  contained  in  the  region  enclosed  by  u  +  v  =  q  , 

v  =  0,  and  rQ ,  where  is  generated  by  (6)  with  =  0  and  ft  >  0. 

Up  to  this  point  the  solution  has  presented  no  particular  difficulties, 

but  equations  (4)  can  only  be  written,  and  hence  the  "solutions"  (6)  and  (7),  if 

u  +  v  -  a  >  0.  Using  the  results  in  appendix  4,  it  is  seen  that  at  the  front, 
2  2  2  2  2 

u+v-a=u  .-a  .  =  0.  This  means  that  the  detonation  front  is  both  a 

+  _  c-j  c-j 

C  and  C  characteristic  (the  equations  are  parabolic).  In  the  theory  of 
characteristics  it  is  shown  that  the  method  can  be  applied  if  the  boundary 
along  which  data  are  given  is  not  a  characteristic,  but  in  this  case  it  is.  The 
way  this  difficulty  was  handled  is  discussed  in  appendix  6. 
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B.  Numerical  Method 


The  equations  to  be  solved  are  (see  appendix  6) 


+  2TT 

r  :  -  =  6  =  R2  cos 

c-j 


()3+  +  0*)  -  R3  cos  [(R2/R3)/3+  +  ipl  ] , 


2V 

A  . 

c-3 


=  X  =  R2  sin  (j3+  +  4k,)  -  R3  sin  (R2/R3)/3  '  +  i //* 


,+  .  .+ 


along  which 


dy/dx  =  -l/tanf  -+  ip* 


5  =  R2  cos  (/ 3"  +i pi)  -  R3  cos  (R2/R3)/3  + 


X  =  R2  sin 


in  (j3~  +  1 pZ)  -  R3  sin  ^(R2/R3)/3_  +  , 


along  which 


dy/dx  =  -1/tan  (j^~J  +  ^*)  * 


where 


0  <  /3+  <  71-d  -  ju)/2ju,  0  >  j3“  >  —  7T ( 1  -  ju)/2ju>  generally, 

+  -  <3) 
and  0  >  y!/;,;  >  7r(l  -  p)/2JuJ  0  <  xfj *  <  7r(l  -  ju)/jn  . 

Specifically  the  limits  on  |3+  and  jB~  are  determined  by  the  T+  curve  with 

du,  -  0  and  X  =  0  as  well  as  by  the  circle  ./ ^ ,  as  given  in  detail 

c-j 

in  appendix  6. 

The  nature  of  the  above  relations  demands  a  numerical  solution.  The 
procedure  will  be  to  assign  values  of  1//^  and  ip ^  which  will  be  used  to  determine 
/3+  and  /3~  and  hence  U  and  V,  then  to  determine  X  and  Y  by  making  numerical 
integration  of  the  compatibility  conditions. 

To  that  end  several  general  relations  will  be  derived  below:  Suppose 
that  one  is  given  1 p+  and  xpZ,  and  is  required  to  find  and  /3+  at  the  intersection 
of  r+  and  T".  Equations  (1)  and  (2)  are  used  to  give 

+  +  /Ro  +  +\ 

R2  cos  (/3  +4/J  -  R3  cosf  /3  +  (//.J  =  R2  cos  (j3  +  ip J 


"  R3  cos(tqP~  +K 
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and 

R2  sinO+  +  (//J)  -  R3  sin  [(R2/R3)0+  +  0*]  =  R2  sin(j3“  +  i^l) 

-  R3  sin  |(R2/R3)/3  +  0,;. 

Eq.  (4)  may  be  rewritten  as 
R 


(4b) 


cos(j3+  +  0*.)  -  cos(/3  +  0.,.) 


r3  < 

[  cos 

J  j 

l  L 

cos 


(R2/R3)^+  +  0* 


(R2/R3)j3'  +  ipl 


(5a) 


and 

R0 


sin(j3+  +  ipt.)  -  sin(/3  +  0;|;) 


:  R3|sin  [(R2/R3)/3  +  0*J 

-  sin  |^(R2/R3)/3  +  0}., 


(5b) 


The  following  relations  are  easily  derived  from  trigonometry: 
cos  * p  -  cos  0  =  2  sin  ^  (0  +  0)  sin  (0  -  0)  , 

sin  0  -  sin  0  =  2  cos  ^  (0  +  0)  sin  (0  -  0)  . 

Using  (6)  in  (5): 

2R„  sin  i  (/ 3+  +  /3  +  0*,  +  0;.; )  sin  (8+  -  0  +  0*  -  0*) 


(6) 


(7) 


=  2R3  sin  tj- 


,+ 


+ 


(R2/R3)(j3  +  j3  )  +  0,.,  +  0.;. 


sm  2 


+ 


+ 


(R2/R3)(]3  -0)+0*-0*J 


and 


2R2  cos  ( P+  +  /3  +  0l  +  0l)  sin  (/3+  -  (3  +  0*1  -  0;|, ) 


(8) 


2R3  cos  ^ 


(R2/R3)(/3+  +  0  )  +  0*  +  0l 


.  1 
sm  2 


(R2/R3)(0+  -  0~)  +0* 


-  0: 


Dividing  equation  (7)  by  (8): 

tan  (0+  +  0  +0t  +  0;|<)  =  tan (R2/R3)(0+  +  j3  )  +  0^  +  0s!< 


(9) 


From  equation  (9): 

+ 


+ 


+ 


0  +  0  +0>;5  +  05|<-  (R2/R3)(]3  +  (3  )+0s!. +  0.,, 
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(1  -  R2/R3)(j3  +/3")  =  0  .  (1 

Recalling  the  definitions  of  R2  and  R3  from  equation  (53),  appendix  6: 

1  -  Vr3  - 1  -  =  t^h  * 0  • 

hence,  (10)  yields 

(3+  =  -/T  .  (1 

Hence,  the  general  statement:  At  the  intersection  of  any  pair  of  F+  and  £ 
characteristics  one  need  only  determine  one  of  /3  and  /3  since  ~  -ft  at 
such  intersection  points. 

Now  putting  (11)  in  (7)  and  (8)  yields: 


2R2  sin  +  sin  4  (2^+  +  ^ -:<)  '  2R3  sin  £(**  + 


X  sin  -| 


2R2  +  +  -\ 


2R2  cos  +  0  *)  sin  |  ( 2/3+  +  ^  =  2R3  cos  ^ 


.  1 
X  sin  — 


from  which 


sin  ft+  +  4(^4  _ 


2R2  +  + 

—  ft  +  0*  -  ip* 

3 


Rq  r  +  1/4- 

sinLRj^  +  si0*  '  ^ 


Equation  (12)  is  then  solved  for  j3  .  As  indicated  in  appendix  6,  along  the 

_|_  _  4* 

boundary  X  =  0,  ip>{,  -  Using  this  relation  in  (12),  the  value  of  ft  along 

X  =  0  is  given  by  the  solution  of 


^/3+  +  4j  ^  °> 


which  is  seen  to  agree  with  setting  v  =  0  in  (1). 

+  -1- 

The  other  boundary  where  j3  is  in  question  is  the  Fq  curve  along  which 
ip+,,  -  0,  The  value  of  jB+  along  r*  is  then  found  from  the  solution  of 
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/-+  1 

_^3 

R2  r>+ 

1  .- 

v3  "  2  ^ 

r2 

sin 

.R3@ 

=  0  .  (14) 

From  equation  (17),  appendix  6,  the  sound  speed  A  may  be  determined  from 


Therefore,  given  a  pair  of  r  curves,  by  specifying  <//t  and  0;,  one  determines 
/3+  from  (13)  or  (14).  With  j3+,  one  can  use  (1)  or  (2)  to  determine  U  and  V 
and  (15)  to  determine  A  so  that  the  problem  is  completely  solved  in  the  U,  V 
plane  since  knowing  A  one  can  find  the  pressure  from  equation  (19),  appendix 
1,  and  the  density  from  equation  (18),  appendix  1,  etc.  The  problem  now  is 
to  determine  X  and  Y  at  each  |3+  in  the  U,  V  plane.  The  known  values  of  X,  Y 
are  (see  Fig.  1 ): 


Y  =  0  along  X  =  0, 

X  =  0,  Y  =  b  along  rj, 

X  =  Y  =  oo  along  U2  +  V2  =  A2  ,/|li2, 

^  J 

X  =  0,  0  <  Y  <  b  at  V  =  0,  U  =  A  . , 

—  —  J 

_l_  2  2  2  2 

0  <  X  <  °o,  b<Y<ooat  the  intersection  of  Fq  and  U  +  V  =  Ac_^/q 

-j-  “h 

The  approach  is  to  assign  a  value  of  let  0i;.  (i)  =  0:  ^ -'(2 )  •  •  •  *  out 

to 


+ 


-7r(l  -  q)  . 
2ju 


determine  values  of  \p~  that  give  the  r“  curves  that  intersect  these  r+  curves 
along  X  =  0  so  that  the  region  is  covered  with  a  mesh  of  r  curves;  determine 
|3+  at  each  intersection  point;  and  then  use  an  average  value  of  the  slopes  in 
(1)  and  (2)  to  obtain  the  X,  Y  point  at  each  intersection. 

At  this  point  it  is  necessary  to  convert  the  equations  for  use  in  the 
numerical  solution.  To  do  this  the  characteristic  curves  are  numbered  with 
the  T+  curves  designated  by  k  =  1,  2, .  .  .  ,  and  the  r~  curves  designated  by 
j  =  l,  2, .  .  .  (see  Fig.  la).  To  do  this,  let 


and 


(16) 


Fig.  1.  The  problem, 
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Since  (3  =  -fi+  at  the  intersection  points,  one  need  only  compute  and  so  let 

fi+  -  Bjj  k  then  U  =  Ik  k,  V  =  Vj>  k,  and  A  =  Aj>  k  .  (17) 

For  convenience  let 


R4  =  l  _  M  =  Rl/R3  ' 

R5  -  R2/R3, 

R6  ~  R3/R2  * 

Then  rewriting  (12),  (13),  (14),  (15),  (1),  and  (2): 


sinh,k+i<wk-  wj>]  -Re sin 


in  [R5Bj,  k  +  i (W 


k  -  W.  )J  -  0  for  a  general 

Bi,k- 


Sln  [Bj,  k  +  Wk  ]  -  R6  Sin  K  Bj,  k  +  Wk]  =  0  for  a  Bj,  k  alone  \  k  =  °-  (1  9> 


sinfe.  ,  -  iw.l  -  R„  sin[R(-B 


J—J  .-i  o  v  v  . 

L  J-k  2  3 


5  j,k  2  'j 


W,  =  0  for  a  B.  ,  along  W1  =  0, 


'j,  k  “  Ac-j  [R2  cos(Bj,  k  +  Wk)  -  R3  cos(R5B.j,  k  +  W 
'j,  k  =  Ac-j  [R2  sin(Bj,  k  +  O  -  R3  sin(R5Bj,  k  +  WD 


^-KZ-^k)2]  (V)' 


Then  one  assigns  values  to  and  for  each  k,  j,  and  calculates  B^  ^  from 

(19)  using  the  Newton- Raphson  Method  (see  appendix  8).  Equations  (20)  are 

then  used  to  determine  U.  ,  ,  V.  ,  ,  and  A.  ,  at  each  B.  ,  . 

J,  kJ  j,  k*  j,  k  j,  k 

In  order  to  determine  x  =  X.  ,  and  y  =  Y.  define 


j+~2>  k 


=  tan--^(B.  k  +  Bj+1_k)  +  l  (W-+W-+1) 


tan-F<Bi,k  +  Bj,k+l>+2  K  +  Wk+1) 
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+ 


and  calculate  F  n  and  F  from  the  given  W,,  W.  and  the  calculated 

i4k  J 

B.  ,  .  The  two  compatibility  conditions  from  (1)  and  (2)  may  now  be  written 
b  k 
as 


-V--Yi,k)ri+i  k  =  W-xi.k(alongr)- 

J  2* 

-  <Yj,  k+l  -  Yj,  k>  Fi  k+l  -  Xj,  k+l  -  Xj,  k  <alo"8  r'>- 


(23) 

(24) 


To  carry  out  the  computation  using  the  known  values  of  x,  y  the  procedure  is 
to  first  obtain  the  x  value  at  k  =  2  using  (24),  then  obtain  x,  y  values  along 
for  each  "j"  using  both  (23)  and  (24).  The  process  is  then  repeated  for  each 
"k."  Referring  to  the  diagrams  in  Fig.  1,  Fig.  lb  represents  the  situation 
along  the  U  axis  and  Fig.  lc  represents  the  situation  in  the  U,  V  plane 
generally.  The  computations  in  these  two  cases  are  the  ones  used  and  will 
be  carried  out  here  as  an  illustration: 

To  determine  the  x  value  along  V  =  0  for  a  given  "k"  value,  assume 

X.  ,  ,  Y.  ,  ,  Y.  ,  ,  ,  =  YMIN  are  known;  then  from  (24), 

1,  k  l,  k  l,  k+l 


Xj,  k+l  =  Xj,  k  "  (Yj,  k+l  '  Yj,  k}  F 


b  k+i 


(25) 


determine  subsequent  values  of  k+1  and  Y^  k+1>  X  and  Y  are  known 
;  j+1,  k;  and  j,  k+l  and  from  (23)  and  (24): 


at 


X3+i,k-xJ,k-n  +  Fj+1;k+|V,k-Fj^_ 

Yj+1,  k+l  =  P  ,  -  F 


"-i  1,4-1 


and 

Xj+1,  k+l  =  Xj+1,  k  ‘  (Yj+l,  k+l  "  Yj+1,  k}  Fj+lj  k+l  * 


(27) 


C.  Example  Problem 

As  an  application  of  this  method,  a  problem  was  solved  for  a  typical 
explosive  utilizing  a  CDC  3600  high  speed  computer.  The  data  used  was  for 
PBX  9404  (see  ref.  4): 
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Chapman- Jouguet  pressure  -  P  .  -  0.39  megabar, 

^  “  3 

Detonation  speed  =  D  =  0.88  cm/qsec, 

3 

Density  of  unburnt  explosive  =  Pq  =  1.84  g/  cm  . 

Using  the  relations  in  appendix  4  the  effective  y  was  calculated  to  be: 

Specific  heat  ratio  =  y  -  2.6535795. 

The  initial  half-thickness  of  the  H.  E.  was  taken  as  0.5. 

The  idea  was  to  assign  a  value  of  zero  to  and  W^,  and  a  fixed  value 
AW  to  be  added  to  W~  to  get  W7+1  and  subtracted  from  to  get  W^+1.  The 
grid  thus  generated  in  the  u,  v  plane  was  to  be  used  to  solve  the  problem. 

While  this  gave  a  rather  finely  spaced  grid  in  the  u,  v  plane  with  AW  =  1°  (or 
0.01745  radian),  the  resulting  grid  in  the  x,  y  plane  had  a  relatively  large 
space  between  the  detonation  front  and  the  first  C+  curve.  When  smaller 
values  of  AW  were  tried,  the  gap  was  seen  to  close  very  slowly  and  it  became 
obvious  that  variable  grid  spaces  would  be  required.  Without  any  attempt  to 

optimize  the  system  of  spacing,  30  C  curves  were  inserted  between  the 

+  -  7 

original  C  lines  and  the  front  ranging  in  spacing  from  AW,  =5X10  X  AW 
o  1 

to  AW„n  =  4  X  10  XAW.  From  there  on  the  spaces  were  equal  with  AW, 

=  2  X  10  1  XAW.  The  solution  was  obtained  with  AW  =  1°  so  the  majority  of 

the  spacing  was  with  0.2°  intervals.  Figures  2  and  3  show  some  representative 

r  and  C  curves  respectively. 

In  considering  the  accuracy  of  the  method  used,  the  C  curve  nearest 
the  front  is  the  most  important.  The  initial  point  on  this  curve  along  y  =  0  is 
determined  by  extending  a  straight  line  from  the  point  x  =  0,  y  =  0.5.  Referr¬ 
ing  to  Fig.  3  it  seems  obvious  that  a  curve  drawn  from  the  point  (0.  ,  0.5)  to 
the  x  axis  without  crossing  the  first  C  curve  cannot  differ  much  from  a 
straight  line.  By  taking  smaller  value  for  AW^,  the  first  point  on  this  C+ 
curve  can  be  placed  as  near  to  the  front  as  desired  and  hence  the  x,  y  values 
can  be  made  as  accurate  as  desired. 

As  a  check,  AW  was  set  to  0.5°  giving  AW ^  =  (2.5  X  10  ^)°  ,  etc.,  and 
the  resulting  values  of  x  and  y  were  the  same  to  four  significant  figures  near 
y  =  0  and  to  less  than  1%  everywhere. 

Figure  4  shows  the  pressure  in  megabars  as  a  function  of  the  distance 
from  the  front  expressed  in  units  of  the  original  thickness  for  various  values 
of  y.  Table  1,  appendix  11,  lists  output  from  which  this  figure  was  plotted. 
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Figures  5,  6,  and  7  show  pressure  in  megabars,  velocity  in  the  x 
direction  relative  to  the  front,  and  velocity  in  the  y  direction  as  a  function  of 
the  distance  from  the  x  axis  for  various  distances  from  the  front. 

Figure  8  shows  constant  relative  volume  contours  in  the  x-y  plane. 

D.  The  Program  Card  Decks 

The  code  name  given  this  program  is  2D-STEDET.  One  basic  program 
was  written  in  FORTRAN  for  the  CDC  3600  high  speed  computer.  There  are 
two  versions  available  and  they  differ  only  in  output  and  COMMON  statements. 
Both  binary  and  FORTRAN  decks  are  available  along  with  copies  of  the 
compilations.  If  the  decks  are  desired,  they  may  be  obtained  by  contacting 
Mark  Wilkins. 

No  attempt  was  made  to  optimize  the  program  to  achieve  rapid  calcula¬ 
tions  or  uniform  distribution  of  grid  points  in  the  x,  y  plane.  This,  of  course, 
could  be  done.  Anyone  desiring  to  use  the  program  will  find  the  following 
comments  of  interest. 

2D-STEDET 

The  program  allows  for  a  7  >  2.6535795  and  gives  as  output  along  each 
j  and  k  curve: 

P  =  pressure, 

VOL  =  relative  volume, 

E  =  internal  energy  in  units  of  original  volume, 

U  =  velocity  in  the  x  direction  relative  to  the  front, 

V  =  velocity  in  the  y  direction, 

B  =  parameter  on  r+  and  r  , 

WPOS  =  generating  parameter  for  F+, 

WNEG  =  generating  parameter  for  F  , 

X  =  distance  from  the  front, 

Y  =  distance  from  the  centerline  of  H.  E. , 

ULB  =  maximum  value  of  B  on  each  k  line. 

In  addition  there  are  linear  interpolation  routines  which  were  put  in  to  get 
information  needed  to  show  the  variations  along  constant  X,  Y,  and  VOL 
lines.  The  first  routine  prints  out  P  and  X  for  Y  =  .  1,  .2,  .3,  .4.  The 
second  prints  out  U,  V,  P,  D-U  =  UPRI,  and  Y  for  X  =  .  25,  .  50,  .  75,  1.0. 

The  third  prints  out  X  and  Y  for  VOL  =  .  8,  .9,  1.0,  1.1,  1.2,  1.3,  1.4,  1.5, 
1.6,  1.7,  1.8. 
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Fig.  6.  Velocity  normal  to  the  front  as  a  function  of  distance  from  the 
centerline  of  the  H.  E.  for  various  distances  from  the  front  in  units  of  the 


initial  slab  thickness,  for  an  ideal  gas  equation  of  state  in  plane  geometry 


With  AW  =  DELWS  =  0.01745  this  takes  about  20  minutes  to  run  to 
completion. 

2D-STEDET  W/OOP 

The  program  was  set  up  to  run  with  7  >  2.65357  95  and  AW  =  0.008725. 

If  AW  =  0.01745,  a  smaller  7  would  be  possible  (see  below).  This  program 
differs  from  2D-STEDET  in  that  the  first  set  of  output  (P,  VOL,  etc.  )  is 
printed  out  only  at  WPOS  =  0.01745  and  WPOS  =  0.1745  for  each  j,  and  in 
addition  to  the  interpolation,  prints  out  P,  D-U,  and  X  for  Y  =  0. 

With  AW  =  0.008725  this  takes  about  6  minutes  to  run  to  completion. 

General  Comments 

In  either  program  if  sense  switch  1  (SSI)  is  depressed,  the  problem 
will  come  off  after  completion  of  calculations  along  a  k-line.  The  on-line 
printer  will  give  the  k-line  number  just  calculated  and  the  maximum  number 
of  k-lines,  and  will  indicate  "PROBLEM  FINISHED"  after  the  output  data  has 
been  written  on  the  output  tape.  Everything  calculated  to  that  point  is  printed 
out  but  the  problem  must  be  resubmitted  if  the  complete  calculation  is  desired. 

The  size  of  the  region  in  the  U,  V  plane  is  dependent  on  7  but  the  grid 
spacing  is  constant  or  at  least  dependent  on  an  unrelated  input  quantity 
DELWS.  As  a  result,  a  change  in  7  results  in  a  change  in  the  number  of  grid 
lines.  The  number  of  j-lines  is  determined  as  follows: 

tt(-  -  l)  2-1 

Vax  -  0.2* DELWS  +  2 9>  M  =  7+~T  * 

In  2D-STEDET  j  should  not  exceed  600  and  in  2D-STEDET  W/OOP  j 

max  max 

should  not  exceed  1200. 

If  the  program  were  to  be  used  often,  it  would  be  advisable  to  establish 

a  fixed  number  of  grid  lines  with  a  grid  spacing  dependent  only  on  7.  As  it  is, 

if  a  problem  is  to  be  run  where  j  is  too  large,  the  COMMON  statements 

max 

in  the  FORTRAN  deck  will  have  to  be  changed.  Care  must  be  taken  in  this 
case  since  2D-STEDET  W/OOP  currently  uses  about  45,000  words  of  memory 
in  the  CDC  3600. 

To  run  a  problem  with  one  of  the  binary  decks  one  needs  only  to  put  a 
*ID  card  in  front  of  the  deck  followed  by  a  *XEQ  card,  the  binary  deck,  a 
*I)ATA  card,  and  three  input  cards.  The  input  cards  must  be  in  the  F10.7 
format  and  in  the  following  order: 
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Card  1: 

DELWS,  YSTRT  =  half,  thickness  of  H.  E.  ,  XSTRT  =  0. ,  WSTRT  =  0.  , 
YM  =  0. ,  ULA  =  0.001 
Card  2: 

GAMMA  =  y,  SPEED  =  D  =  detonation  velocity,  ROE  =  pQ  =  reference 
density  of  the  H.  E. 

Card  3: 

CY1  =  0.1,  XC1  =  0.25,  CVOL1  =  0.8 

The  three  input  quantities  on  Card  3  are  used  to  activate  the  interpola¬ 
tion  routines.  CY1  gives  P  and  X  for  a  fixed  Y;  XC1  gives  U,  V,  P,  D-U 
=  UPRI,  and  Y  for  a  fixed  X,  CVOL1  gives  X  and  Y  for  a  fixed  VOL.  If  any 
of  these  is  not  desired,  then  the  corresponding  activation  number  (CY1,  etc.) 
should  be  put  on  Card  3  as  a  "0."  since  the  program  is  set  to  bypass  for  that 
value. 

2D-STEDET  W/OOP  PLUS  PLOT 

This  program  is  the  same  as  2D-STEDET  W/OOP  except  that  a  plot 
routine  has  been  added  which  will  plot  the  output  obtained  from  the  interpola- 
tion  routine  plus  every  twentieth  T  and  C  curve.  The  r  curves  start  with 
k  =  1  while  the  C+  curves  start  with  k  =  2.  The  plot  routine  can  be  bypassed 
completely  or  the  r+  and  C+  curves  may  be  left  out  in  addition  to  the  bypass 
already  available  for  the  interpolation  routines. 

To  use  this  program  with  a  binary  check  an  additional  card  is  required 
in  the  110  format: 

Card  4: 

IPLOT  -  IA,  IPLOT1  =  IB 

If  IA  =  0  and  IB  =  0,  no  plot  results. 

If  IA  =  1  and  IB  =  0,  C+  and  r+  are  not  plotted. 

If  IA  =  1  and  IB  =  1,  C+  and  r+  are  plotted  and  the  interpolation  curves  are 

calculated  and  plotted  according  to  the  data  on  Card  3. 

An  example  of  these  plots  is  shown  in  Figs.  2-8. 
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II.  A  PLANE  TWO-DIMENSIONAL  DETONATION  FOR  AN 
EXPLOSIVE  WITH  A  "WILKINS"2  EQUATION  OF  STATE 

Introduction 

In  parts  A  through  D  below  is  presented  the  solution  of  a  plane  two- 
dimensional  steady- state  detonation  of  an  explosive  with  a  "Wilkins-type" 
equation  of  state  employing  the  method  of  characteristics.  Part  A  gives  an 
outline  of  the  derivation  of  the  basic  equations;  part  B  gives  a  description  of 
the  numerical  method  used;  part  C  describes  the  example  problem  solved; 
and  part  D  gives  a  description  of  the  program  card  decks  available. 


A.  Calculation  of  the  Two-Dimensional  Steady-State  Detonation 
for  an  Explosive  Having  a  "Wilkins"  Equation  of  State 


The  equation  of  state  used  here  is  that  given  in  ref.  4.  The  basic  equa¬ 
tions  of  motion  are  the  same  as  those  in  section  I,  part  A,  Eq.  (1).  The 
difference  comes  in  Eq.  (2),  section  I,  part  A,  which  is  replaced  by  (see 
appendix  9): 


u 


2 


+  v 


2 

—  dp  =  const. 
P 


(1) 


By 


definition  a 


2  _  9p 
=  ^P 


,  and  from  the  Wilkins  equation  of  state: 
r? 


p0 


Ag^  +  BR  vVRV  + 


(1  +  w)C 


V 


u 


(2) 


When  (2)  is  substituted  in  (1): 


2  ,  2.2 

u  +  v  +  — 

p0 


AQ 


(Q  -  1)VQ_1 


B(v  +  i) 


-RV 


(1  +  u)C 


uV 


w 


const.  (3) 


The  constant  in  (3)  can  be  evaluated  at  the  detonation  front  where  v  =  0, 

u  =  a  .,  V  =  V  ..  This  implies  that  through  (3)  the  relative  volume  (V)  is 
J  ^  “  J  2  2 

a  function  of  u  and  v  and,  since  a  =  a  (V),  that  a  is  a  function  of  u  and  v. 
Hence,  the  hodograph  transformation  can  be  made  for  this  set  of  equations 
just  as  in  section  I,  part  A,  and  in  fact  all  the  results  through  appendix  5  are 
unaffected.  Hence,  we  have: 
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r  :  v  =  tan  d/  u  , 
a  *  a 


+ 

C  :  y. 


tan  ip 


r":  -  tan  \p  , 

C":  Jo  =^-+x«  • 
p  tan  ip  p 

The  problem  from  here  is  different.  The  solution  for  T  and  T  cannot 
be  found  in  closed  form  as  in  section  I,  but  a  numerical  integration  must  be 
performed  to  obtain  the  values  of  ip+  and  ijj  ,  the  slopes  of  the  tangents  of  T 
and  r"  respectively  (see  appendix  10).  Having  determined  these  parameters, 
the  equations  for  C+  and  C_  can  be  used  to  determine  x  and  y.  Otherwise, 
the  solution  is  very  similar  to  the  ideal  gas  equation  of  state  problem. 

B.  Numerical  Method 

Referring  again  to  Fig.  1  c,  the  equations  (4)  of  part  A,  section  I,  can 
be  put  in  difference  form  as: 

Yj+1,  k+1  ‘  Yj,  k+1  =  TF2  (Xj+l,  k+1  "  Xj,  k+l)} 


Yj+1,  k+1  “  Yj+1,  k  =  TF1  (Xj+l,  k+1  "  Xj+1,  k)j 


which,  when  solved  for  Xj+1  k+1>  Yj+1  k+1’ 

v  Xi,k+1  -  xi+l,  k  +  TF2Yi.  k+1  -TF1Y1+l,k 

X  j+l,  k+1  TF2  -  TF1 

Xj+1,  k+1  =  Xj+1,  k  +  TF1*Yj+l,  k  "  Yj+1,  k+1  ^ 


where 


TF1  =  tan 


2  C^j+l,  k+1  +  ^j+l,  k)  * 

2(^1+!.  k+1  +  k+l)  * 


TF2  =  tan^0j+1>k+1  +^k+i;  • 

The  equations  involving  u  and  v  are  solved  in  appendix  10  and  yield  the 
following: 


-29- 


+ 


-r= 

Jv, 

* 


+ 


F(v)dV  +  xp , 


(4) 


V 


=  f  F(v)d  V +<//", 
JV„. 


where 


1  da 


F(v) 


_  2  dV 


a 

V 


ga 


(5) 


and  g,  a,  a^  are  all  known  functions  of  V. 

In  order  to  begin  the  solution  it  is  necessary  to  integrate  (4)  which  is  an 

improper  integral.  Using  the  knowledge  gained  in  section  I,  it  is  seen  that 

the  r+  curve  corresponding  to  the  uppermost  point  on  the  detonation  front 

(r*)  has  xp+,.  =  0  and  V*  =  V  ..  At  V  =  V  ,,  g  =  0  and  hence  care  must  be 

taken  with  (4).  The  evaluation  procedure  consisted  of  three  steps:  (1)  the 

integral  was  shown  to  exist  for  V*  =  V  .in  appendix  10,  (2)  an  estimate  of 

c-j 


ip+  for  a  value  of  V  slightly  greater  than  V  .  (0.7263  for  PBX  9404  and 

0.731575  for  LX-04-01)  was  obtained  by  fixing  the  upper  limit  and  varying  the 

lower  limit  toward  V  .  as  far  as  possible  and  then  extrapolating  to  V 

(3)  and  finally,  the  value  of  this  first  xp  was  arbitrarily  varied  to  give  a 

positive  value  of  the  y- velocity  (v)  which  is  required  physically.  The  latter 

was  possible  because  the  other  quantities  calculated  were  not  very  sensitive 

to  small  changes  in  the  small  first  value  of  xp+.  The  integration  was  performed 

5 

using  the  subroutine  ROMBRG. 

With  this  first  value  of  0+  calculated,  arbitrary  values  of  V  were 
assigned  up  to  V  =  10  and  ip+  was  calculated  for  each  V  from  (4).  Equation 
(15),  appendix  10  ,  was  then  used  to  calculate  ip  for  each  V.  This  gives 
everything  on  since  x  and  y  are  constant  along  this  curve  as  shown  in 
section  I. 

The  value  of  V*  for  the  next  T+  curve  is  found  using  equation  (18), 
appendix  10,  and  calculating  ipt,,  and  ip from  equations  (16)  and  (17),  appendix 
10.  Equation  (19),  appendix  10,  is  used  to  calculate  V  for  succeeding  points 
on  this  r  curve,  and  then  equation  (15)  is  used  to  calculate  xp  .  Equations 
(13),  appendix  10,  was  used  to  calculate  u  and  v  and  equations  (14)  and  (16) 
in  UCRL-7797  (reproduced  below)  were  used  to  calculate  p  and  E.  With 
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values  of  ip+,  ip  ,  a,  g  on  two  k  lines,  X  and  Y  can  be  calculated  from  (2) 
above.  The  process  is  then  repeated  for  the  next  k  line. 

C 

nr\  A  ,  dp'RV  ,  s 
p(v)  =  -3  +  BC  +  ^7  ■ 

C  ,^_i>  e-BV_ 
w  uV^'1  w 

where  C,  A,  B,  R,  C  ,  w  are  constants  defined  in  part  D  below. 

S 


C.  Example  Problem 

As  an  application  the  same  problem  was  solved  as  that  discussed  in 

section  I.  The  results  are  shown  in  Figs.  9-15. 

■f* 

Figure  9  shows  some  r  curves. 

Figure  10  shows  some  C+  curves. 

Figure  11  shows  pressure  in  megabars  as  a  function  of  the  distance 
from  the  front  expressed  in  units  of  the  original  thickness  for  various  values 
of  y.  Table  2,  appendix  11,  lists  output  from  which  this  figure  was  plotted. 

Figures  12,  13,  14  show  pressure  in  megabars,  velocity  in  the  x 
direction  relative  to  the  front,  and  velocity  in  the  y  direction  as  a  function  of 
the  distance  from  the  x  axis  for  various  distances  from  the  front. 

Figure  15  shows  constant  relative  volume  contours  in  the  x,  y  plane. 


D.  The  Program  Card  Decks 

The  code  name  given  to  this  program  is  2D-STEDET  —  W/WEOS.  One 
program  was  written  in  FORTRAN  for  the  CDC  3600  high  speed  computer. 
Both  binary  and  FORTRAN  decks  are  available  along  with  copies  of  the 
compilations.  If  the  decks  are  desired,  they  may  be  obtained  by  contacting 
Mark  Wilkins.  As  in  the  ideal  gas  equation  of  state  case,  no  attempt  was 
made  to  "tidy  up"  the  program. 

The  output  available  consists  of: 

SCJ  =  Chapman- Jouguet  sound  speed, 

D  =  constant  in  equation  (3),  appendix  10, 

V  =  relative  volume, 

XVEL  =  velocity  in  the  x  direction  relative  to  the  front, 

YVEL  =  velocity  in  the  y  direction. 
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Fig.  11.  Pressure  profiles  at  various  distances  from  the  centerline  of 
the  H.  E.  as  a  function  of  distance  from  the  front  in  units  of  the  initial  slab 


thickness,  for  the  Wilkins  equation  of  state  in  plane  geometry. 


Y-AXIS 


Fig.  13.  Velocity  normal  to  the  front  as  a  function  of  distance  from  the 
centerline  of  the  H.  E.  for  various  distances  from  the  front  in  units  of  the 
initial  slab  thickness,  for  the  Wilkins  equation  of  state  in  plane  geometry. 


YVEL-AXIS 


Y-AXIS 


Fig.  14.  Velocity  parallel  to  the  front  as  a  function  of  distance  from  the 
centerline  of  the  H.  E.  for  various  distances  from  the  front  in  units  of  the 
initial  slab  thickness,  for  the  Wilkins  equation  of  state  in  plane  geometry. 
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P  -  pressure, 

E  =  internal  energy  in  units  of  the  original  volume, 

-f~  ”f“ 

PSIP  =  angle  ip  made  by  r  with  v  =  0, 

PSIN  =  angle  i p~  made  by  r  with  v  =  0, 

X  =  distance  from  front, 

Y  =  distance  from  centerline  of  H.  E. , 

UPRI  =  velocity  in  the  x  direction, 

U  =  XVEL. 

In  addition  there  is  a  linear  interpolation  routine  exactly  like  that  des¬ 
cribed  for  2D-STEDET  W/OOP  PLUS  PLOT. 

To  run  a  problem  with  one  of  the  binary  decks  one  needs  only  to  put  a 

*ID  card  and  a  *XEQ  in  front  of  the  deck,  the  binary  deck,  a  *DATA  card, 

and  five  input  cards.  These  cards  are: 

Card  1:  5F  10.7  format 
C,  Q,  R,  B,  W 

The  five  quantities  are  the  constants  a,  Q,  R,  B,  u  that  appear  as 
equation  (15)  in  UCRL-7797. 

Card  2:  5F  10.7  format 

ROE,  CS,  VCJ,  SPEED,  YSTRT 
These  five  quantities  are: 

Pq  =  ROE  =  reference  density  of  the  H.  E.  , 

CS  =  Cg  =  the  constant  referred  to  at  the  top  of  page  12  in 
UCRL-7797, 

VCJ  =  Chapman- Jouguet  reference  volume, 

SPEED  =  detonation  velocity, 

YSTRT  =  half-thickness  of  the  H.  E. 

Card  3;  3F  10.7  format 

CY1  =  0.1,  XC1  =  0.25,  CV1  =  0.8 

Card  4:  31  10  format 

ISTRT  =  IC,  IPRT1  =  ID,  IPRT2  =  IE 
where: 

IC  =  0  for  PBX  9404, 

IC  =  1  for  LX-04-01, 

ID  =  1  yields  printout  of  XVEL,  YVEL,  P,  E,  V  for  each  j,  k. 


-3  8- 


ID  =  0  yields  no  printout  of  XVEL,  YVEL,  P,  E,  V, 

IE  =  1  yields  printout  of  X,  Y,  V  for  each  j,  k, 

IE  =  0  yields  no  printout  of  X,  Y,  V. 

Card  5:  21  10  format 

IPLOT,  IPLOT1  (See  explanation  for  2D-STEDET  W/OOP  PLUS  PLOT) 
This  program  will  run  a  problem  for  either  PBX  9404  or  LX- 04 -01.  The  data 
required  for  each  is: 


PBX  9404 

LX-04-01 

C  =  -0.004563 

C  =  -0.0008335 

Q  =  4.0 

O 

U 

G> 

R  =  4.0 

R  =  4.0 

B  =  6.572 

B  =  5.943 

W  =  0.35 

W  =  0.4 

ROE  =  1.84 

ROE  =  1.865 

CS  =  0.032 

CS  =  0.029 

VCJ  =  0.7262958 

VCJ  -  0.73156  94 

SPEED  =  0.88 

SPEED  =  0.84  8 

ISTRT  =  0 

ISTRT  =  1 
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III.  A  CYLINDRICAL  TWO-DIMENSIONAL  DETONATION 
FOR  AN  EXPLOSIVE  WITH  A  "WILKINS" 

EQUATION  OF  STATE 

Introduction 

In  parts  A  through  D  below  is  presented  the  solution  of  a  cylindrical 
two-dimensional  steady-state  detonation  of  an  explosive  with  a  "Wilkins" 
equation  of  state  employing  the  solution  to  the  plane  problem  discussed  in 
section  II  to  start  the  solution,  and  then  employing  the  method  of  character¬ 
istics  for  the  cylindrical  problem  from  there  on.  Part  A  gives  an  outline  of 
the  derivation  of  the  basic  equations  governing  the  motion;  part  B  gives  a 
description  of  the  numerical  method  used;  part  C  describes  the  example 
problem  solved;  and  part  D  describes  the  program  card  decks  available. 


A.  Calculation  of  a  Steady-State  Detonation  in  Cylindrical  Coordinates 
for  an  Explosive  with  a  "Wilkins"  Equation  of  State 


In  considering  the  cylindrical  case  it  is  only  necessary  to  note  that 
there  is  no  change  in  the  basic  derivation  given  in  section  II  except  that  the 
conservation  of  mass  equation  is  altered  to: 


9p ,  9 p  ,  9 v  ,  9u_  v 
Fy  U  Fx  ^  Fy  ^  0x  "  ^  y 


(1) 


where  x  =  coordinate  along  the  axis  of  the  cylinder, 

y  =  coordinate  normal  to  the  axis  of  the  cylinder, 
and  it  is  assumed  that  there  is  no  rotation  about  the  axis.  The  other  con¬ 
servation  equations  and  Bernoulli' s  equation  remain  the  same. 

Following  the  same  procedure  used  in  section  II,  the  set  of  equations  to 
be  solved  is: 


/  2  2, 

(a  -  u  ) 


9u 

F3F 


9v 

UVF7 


9u 

UV  ts—  + 

9y 


(a‘ 


v2)p-* 


o 


9y  9u 
9x  Fy 


+  0 


(2) 


2 

cl  V 

Due  to  the  presence  of  the  term - on  the  right  of  (2),  the  hodograph  trans- 

y 

formation  used  in  section  I  will  no  longer  yield  a  set  of  uncoupled  equations 
in  characteristic  velocity  space.  The  derivation  of  the  characteristic  equa¬ 
tions  in  the  x,  y  plane,  however,  follows  the  same  procedure  but  this  time  the 
vector 
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C  - 


r  2 
I  a  v 

L  y 


(3) 


With  this  value  of  C,  a  straightforward  application  of  the  technique  followed 
in  appendix  2  yields: 


C 


.  dy 

1*  dx 


a 


uv 


(4) 


along  which 


du 

+  1 
a 

f  uv  +  a  /u2  +  v2  -  a2%N 

j  dv 

2 

a  v 

dx 

k  u2  -  a2  J 

'  dx 

a  y(u  -  a  ) 

and 


C  •  rJL 
^ 2 •  J~ 


dy 

dx 


,  /T 

uv  +  a  vu 


+  v 


P 


u 


(5) 


along  which 


du 

dx 


+ 


uv 


/ 2  ,  2  2\  , 

-  a  vu  +  v  -  a  \  dv 

)  dx 


P 


2  2 
u  -  a 


2 

a  v 


J3  y(u2  -  a2) 


P 


In  addition,  the  following  relations  from  section  II  obtain: 
2  1 


0 


V 


AO  2  -RV  ^  +  U^s 

+  BRV2e  RV  +  — 3-5 


2  2  ,  2  2 

q  =  u  +  v  =  D - 

P 


0 


V 


AQ  _H-B(v^)e-RV  +  -— 


(Q  -  l)V 


W- 


coV 


10 


and 


~  2  ,  2 
D  =  a  .  +  — 
e-J  P0 


c-j 


-RVc_j  i  (1  +u)Cs 


+ 


-r  rlO 

10V 

c-J 


(6) 


(7) 


(8) 


In  order  to  start  the  solution,  it  is  assumed  that  "near"  the  front  (small 
x)  the  solution  is  approximately  the  same  as  that  for  the  plane  case,  and 
hence,  approximate  values  of  u,  v,  V  and  therefore  p,  a,  E  at  specified 
values  of  y  on  an  x  line  "near"  the  front  can  be  obtained  from  the  solution 
with  2D-STEDET  W/WEOS  as  described  in  section  II.  This  assumption 
appears  reasonable  for  a  sufficiently  small  x,  but  no  attempt  has  been  made 
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to  check  the  assumption  or  even  to  vary  the  value  of  the  "small  x,"  due  to 
lack  of  time.  The  program  was  written  with  the  latter  in  mind,  and  it  could 
easily  be  done. 

In  order  to  solve  equations  (4)  and  (5),  use  is  made  of  the  fact  that 
along  the  x  axis,  y  =  0  and  v  =  0.  Since  the  term  v/y  occurs  in  (2),  one  can¬ 
not  simply  set  y  and  v  to  zero  in  (4)  and  (5).  Taking  the  limit  of  the  first  of 
(2)  as  y  and  v  tend  to  zero  yields: 


(a2- 


2. 
u  ) 


3u 

35 


-  0 


9u 

3v 


+ 


2  9v 

1  *y 


=  -  a 


2  9v 
■gy  j 


(9) 


since  v/y  is,  in  fact,  (v  -  0)/(y  -  0),  and  as  y  tends  to  zero  is,  by  definition. 
With  (9)  and  the  second  of  (2)  in  place  of  (2),  the  characteristic  equations 
along  y  =  0  are: 


(10) 


along  which 


du  ,  / 2  a  dv 

dx  2  2  dx 

a  u  -  a 


a 


=  0  (along  y  =  0)  , 


and 

dy  _  {2  a 

z’  dx  0  2  2 

/3  u  -  a 


along  which 


du 

dx 


0 


/2  a  dv 

~2  2  dx 

u  -  a 


/3 


=  0  (along  y  -  0)  . 


(11) 


Equations  (4)  through  (8),  (10),  and  (11)  along  with  the  equations  for  p 
and  E  from  part  B,  section  II,  constitute  the  set  of  equations  required  for  the 
solution  of  the  cylindrical  steady- state  detonation  problem. 


B.  Numerical  Method 

The  numerical  procedure  is  rather  straightforward.  The  idea  is  to 
replace  the  characteristic  equations  of  part  A,  section  III,  with  difference 
equations,  to  use  values  for  all  the  variables  along  a  selected  x  =  constant 
line  determined  from  the  program  described  in  section  II  for  generating 
values  on  an  adjacent  curve  in  the  "cylindrical  region,"  and  to  repeat  the 
generating  process  until  the  problem  is  finished.  To  that  end  define: 
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T 


T 


+ 

T 


+ 

T 


UV 


-  a/u2  +  v2  -  a2 


1 - “T 

u  -  a 


-/2  a 

“2 - 2 

u  -  a 


.  /  2  .  2  2 
uv  +  avu  +  v  -  a 

- 2 - 5 - 

u  -  a 

/2  a 

“2 - 2 

u  -  a 


in  general, 
on  y  -  0, 


in  general, 

on  y  =  0, 


(1) 


(2) 


e  = 


2 

a  v 

"2 - 2 

u  -  a 


0=0 


in  general. 


on  y  =  0, 


(3) 


AG  =  -  Gj  for  i,  j  =  1,  2,  or  3  and  G  =  any  of  the  variables  x,  y,  etc. 

Then  the  characteristic  equations  (4),  (5),  (10),  and  (11)  in  part  A,  section 


III,  become: 


C  :  Ay  =  t ”  Ax  along  which  Au  +  r+Av  =  ~y~~ 
C+:  Ay  =  r+Ax  along  which  Au  +  r~Av  =  —  ~  , 


(4) 


If  Fig.  16  represents  the  general  situation  in  which  all  values  of  the 
variables  concerned  are  known  at  points  1  and  3,  then  the  equations  (4)  can 
be  written  as  follows: 
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on  C  :  y2  "  ^3  =  \  2  2  <X2  "  X3*  1 

u2  "  U3  +  2  2 2  ^V2  "  V3)  =  2  6 2  *X2  '  X3)  * 

On  C  ;  ^2  ~  ^  \  ~2  ^2  x]^  > 


v 2  "  ui  +  ^  (v2  ~  ^ ~  \  6i  ^x2 


x1)  . 


In  these  equations, 

2^  =  t g  +  with  p  =  +  or 


(5) 


(6) 

(7) 

(3) 


(9) 


S 


P 

2 


with  p  =  +  or 


(10) 


6 


1 


(11) 


(12) 


Equations  (5)- (8)  constitute  a  set  of  equations  to  be  solved  for  the  values  at  2, 
ignoring  for  the  moment  the  problem  of  evaluating  r  ,  t  and  6  at  2.  Solving 
(5)  and  (7)  yields: 


2(\-  -  v  )  -1-  ^  x  -  x 

‘v  3  .  ]_ <  —  ^  -'q  —9  a3 


.xf 


(13) 


anti 


>  2  y  1 4”  2  “1  ^x2  ~  xi^’ 


•  (14) 


where 


\  —  v"^"  V  *" 

1  “1  "  “2  ' 


(15) 
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and  the  value  of  x 2  from  (13)  is  used  in  (14).  Solving  (6)  and  (8)  yields 

2(u3  -  Ul)  "  V1  +  S2  v3  +  62(x2  "  x3}  "  51(X2  "  xl} 

V2  =  ^ 


(16) 


and 


tu  =  u,  - 


(V2  - 


V  +  2  61(X2  - 


Xl} 


where 


(17) 


(18) 


and  the  values  of  x2  and  from  (13)  and  (14)  are  used  in  (16)  and  (17)  and 
the  value  of  v„  from  (16)  is  used  in  (17). 

In  order  to  start  the  solution,  the  values  of  t  ,  r  ,  and  6  at  2  are  set 
to  zero.  The  value  of  the  quantities  in  (9)  through  (12)  are  multiplied  by  2 
and  a  first  estimate  for  x2,  y2,  v2,  and  u2  is  determined  from  (13),  (14), 

(16),  and  (17).  To  get  a  second  estimate  it  is  necessary  to  determine  a  value 
for  a  at  2.  To  determine  a2,  it  is  necessary  to  determine  V  at  2.  To  that 
end  equation  (7),  part  A,  section  III,  is  used  as  follows. 

Assume  the  values  of  u2  and  v2  as  calculated  in  the  first  estimate  are 
correct,  then  one  must  solve 


<MV)  =  u2  +  v2 


D  + 


AQ -  i  B(V  -  )e‘ 


(Q  -  1)VQ_1 


+ 


wV 


w 


=  0.  (19) 


This  will  be  solved  using  the  Newton-Raphson  method  described  in  appendix 
8.  If  V1  is  the  ith  estimate  of  the  V  which  solves  4*(V)  =  0,  then 


V1  =  V 


i-1 


1  + 


2a2(Vi_1) 


V 


i-1 1 


< 


(20) 


e,  where  e  is  an  arbitrary  small 


Equation  (20)  is  applied  until  lv* 

_  O 

number  set  equal  to  10  in  this  case.  Equation  (6),  part  A,  section  III,  is 
used  to  calculate  a2(V^  ^ ). 

One  is  then  in  a  position  to  calculate  values  of  r  ,  r  ,,  and  9  at  2. 
Having  done  so,  the  entire  process  is  repeated  to  get  new  estimates  of  x,  y, 
u,  v,  and  V  at  2.  It  was  decided  (for  no  particular  reason)  that  when  two 
successive  values  of  u ^  differed  by  10  or  less,  then  the  values  of  all  vari¬ 
ables  at  2  were  sufficiently  accurate. 
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If  the  point  2  is  on  y  =  0,  then  from  (5)  and  (6): 


x 


2 


^3 


2 


2 


(21) 


and 


Ur 


,  1 
u3  +  2 


22  V 


3+2  ^2^X2  x3^ 


(22) 


where  the  values  of  r~ ,  t+ ,  6  for  2  are  determined  using  the  expressions  in 
(1),  (2),  and  (3)  on  y  =  0.  If  point  1  is  on  y  =  0,  then  there  is  no  change  from 
(13)  through  (18)  except  in  the  calculation  of  r“,  t+,  and  0  at  1  as  indicated 
for  y  =  0. 

In  the  process  just  described  it  is  obvious  that  if  one  has  all  the  vari¬ 
ables  determined  on  some  line  x  -  x  for  N  values  of  y,  then  the  first  set  of 

s 

calculations  will  give  all  the  variables  at  N  -  1  new  points.  This  will,  after 

N  sets  of  calculations,  reduce  the  number  of  new  points  to  1.  However,  if 

the  first  point  on  x  is  used  in  conjunction  with  the  line  y  =  0  through  (21)  and 

(22),  then  a  new  point  is  added  on  the  first  set  of  calculations,  and  on  every 

other  set  of  calculations  thereafter.  The  result  is  that  new  points  are  reduced 

to  1  after  2N  -  1  sets  of  calculations  instead  of  after  N  -  1  with  the  last  point 

being  on  y  =  0.  (In  the  program  written  the  problem  terminates  after  2N  -  2. ) 

Hence,  the  area  covered  by  the  solution  is  determined  by  the  number  of  points 

brought  forward  from  the  initiating  plane  solution  to  the  line  x  =  xg  and  to 

some  extent  upon  the  distribution  of  these  points.  In  the  program  as  written 

no  effort  was  made  to  select  an  "ideal"  distribution  of  points  onx:  x  ,  and 

in  fact,  an  arbitrary  number  of  points  was  selected  which  will  vary  with  the 

equation  of  state  in  an  unknown  manner.  The  distribution  of  these  points  is 

equally  arbitrary.  The  object  of  the  program  was  to  demonstrate  that  the 

method  works.  A  clever  programmer  could  devise  a  routine  to  give  a 

constant  number  of  points  on  x  =  xg  with  an  "ideal"  distribution.  Further, 

the  value  of  xg  was  set  equal  to  0.05  for  no  particular  reason.  This  number 

(x  )  should  be  determined  from  experimentation  with  the  program  in  conjunc- 
s 

tion  with  physical  experiments  and/or  other  computer  programs  such  as 
HEMP  (see  ref.  4). 
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C.  Example  Problem 

As  an  application  of  this  method  the  same  problem  was  solved  as  that 
discussed  in  section  I.  The  results  are  shown  in  Figs.  17-19. 

Figure  17  shows  the  C"  curves  in  the  part  of  the  problem  where  the 
cylindrical  equations  were  solved.  The  blank  space  to  the  left  of  the  line 
x  =  0.05  represents  the  part  of  the  space  where  the  plane  solution  was 
assumed  to  hold.  The  points  calculated  there  were  not  plotted  because  it  was 
felt  they  would  only  confuse  the  figure.  It  will  be  noted  that  the  distance  be¬ 
tween  the  characteristic  curves  near  the  top  of  the  area  becomes  quite  large 
toward  the  end  of  the  solution  area.  This  distance  is  probably  too  great  for 
the  straight  line  approximations  made.  However,  the  addition  of  several  more 
points  would  eliminate  this  objection.  Note  also  that  there  is  a  relatively 
large  gap  between  the  highest  point  on  the  figure  and  the  line  y  -  0.5  which 
later  represents  the  original  boundary  of  the  H.  E.  This  gap  could  be  reduced 
by  the  addition  of  more  curves  in  the  solution  to  the  left  of  x  =  0.05.  The 
curves  are  shown  merely  as  an  illustration  of  the  type  of  information  pro¬ 
duced  by  the  code. 

Figure  18  shows  the  pressure  profile  in  megabars  as  a  function  of  the 
distance  from  the  front  expressed  in  units  of  the  initial  thickness  of  the  H.  E. 
for  various  values  of  y  as  plotted  by  the  plot  routine  directly  on  the  CRT 
(cathode  ray  tube)  by  the  computer. 

Figure  19  shows  the  same  information  as  Fig.  18,  with  the  points  for 
the  plane  solution  (x  <  0.05)  removed,  hand-plotted  on  a  different  scale  for 
clarity.  Table  3,  appendix  11,  lists  output  from  which  this  figure  was 
plotted. 


D.  The  Program  Card  Decks 

The  code  name  given  to  this  program  is  CYL-STEDET.  One  program 
was  written  in  FORTRAN  for  the  CDC  3600  high  speed  computer.  Both 
binary  and  FORTRAN  decks  are  available  along  with  copies  of  the  compila¬ 
tions.  If  decks  are  desired,  they  may  be  obtained  by  contacting  Mark 
Wilkins. 

The  output  quantities  are  the  same  as  those  listed  in  part  D,  section  II, 
with  the  addition  of 

S  =  sound  speed. 
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Fig.  18.  Pressure  profiles  at  various  distances  from  the  centerline  of 
the  H.  E.  as  a  function  of  distance  from  the  front  in  units  of  the  initial 
thickness,  for  the  Wilkins  equation  of  state  in  cylindrical  geometry. 
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The  linear  interpolation  routines  are  the  same  as  in  the  two  plane  cases 
(see  note  below). 

To  run  a  problem  with  the  binary  deck,  follow  the  instructions  given  in 
part  D,  section  II,  with  the  following  change: 

CARD  5:  3110  format 

IPLOT,  IPLOT1,  IPLOT2 

IPLOT  and  IPLOT1  are  unchanged  from  the  plane  cases. 

IPLOT2  =  0  yields  no  plot  of  the  characteristics  in  the  plane  part  of  the 
solution. 

IPLOT  2  =  1  yields  a  plot  of  the  characteristics  in  the  plane  part  of  the 
solution. 

It  is  also  possible  to  solve  a  problem  for  an  ideal  gas  equation  of  state 
for  PBX  9404.  The  following  input  data  is  required  to  do  this: 

C  =  0. 

Q  =  4.0 
R  =  4.0 
B  =  0. 

W  =  1.65358 
ROE  =  1.84 
CS  =  0.1668537 
VCJ  =  0.7262958 
SPEED  =  0.88 
ISTRT  =  2 

When  this  input  was  tried,  the  problem  ran  to  completion,  but  due  to 
the  fixed  selection  of  curves  in  the  plane  part  of  the  solution,  the  maximum 
value  of  x  obtained  was  0.258.  That  was  not  enough  for  a  useful  solution. 

The  important  thing  is  that  this  problem  ran. 

NOTE:  Output  curves  corresponding  to  Figs.  12-15  are  not  displayed  for 
this  program  because  the  area  covered  by  the  solution  obtained  was  too  small 
to  include  more  than  a  few  points  on  x  =  0.25,  0.50,  0.75,  and  1.0  and  almost 
no  points  on  the  constant  volume  contours. 
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APPENDIX  1 

DERIVATION  OF  THE  TWO-DIMENSIONAL,  STEADY,  IRROTATIONAL, 
ISENTROPIC  HYDRODYNAMICS  EQUATIONS 


The  general  equations  governing  continuous  media  are  (ref.  6,  Ch.  2, 


3,  4): 

Conservation  of  Mass 


dp 


w 


+ 


=  0, 


Conservation  of  Momentum 


t 


ij 

;i 


DuK 

Dt  ) 


--  o, 


t^  =  t^. 


Conservation  of  Energy 


P 


Drj 

Dt 


>  0 


Entropy  Principle 


where 

p  5  mass  density, 
u  s  velocity  vector, 

t1^  e  stress  tensor, 
f  i  body  force  vector  per  unit  mass, 
e  =  internal  energy  per  unit  mass, 
p  =  entropy  per  unit  mass, 

A.\  =  covariant  partial  derivative  of  A1, 

-7T  e  material  derivative. 


If  one  considers  a  nonviscous,  nonconducting  medium  in  which  u 
=  [u,  v,  0]  ,  p  =  p(x,  y),  p  =  p(x,  y),  u  =  u(x,  y),  e  =  e  (x,  y),  r)  =  r)(x,  y),  f1  e  0, 
and  t1^  =  -pS1^,  then  these  equations  become: 
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Conservation  of  Mass 


9(pu)  ,  9(pv) 

9x  9y 


0, 


Conservation  of  Momentum 


9p  9u  ,  9u 
cJx  _  P  u  Fx  +  P  v  Fy 


9p  9v  ,  9v 
“F?=pU^+pV^ 


Conservation  of  Energy 


pu 


9x 


lT\  4.  3  (  ,  U‘ 


8 7  (pu)  -  (pv). 


Fx 


where 


Entropy  Principle 


9 n  ,  3p  n 

pu  F^  pV  3y-  °' 


p  £  hydrostatic  pressure, 
x,  y  =  rectangular  coordinates, 

TT2  -  2  .  2 

U  su  +  v  . 

The  equation  of  state  for  an  ideal  gas  is 


where 


P  - 


(  — ) 
\P0J 


,7  /h  "  r?n\ 

exp(— - - )  (see  Ref.  2,  p.  6  ff) 


v.. 

7  s  -£■  =  gas  constant,  or  specific  heat  ratio. 


c  £  specific  heat  at  constant  pressure, 
P 

cv  h  specific  heat  at  constant  volume. 


The  sound  speed  is  given  by 

a2  -  8p 
a  Vp 


so  that  using  (6): 

a2  -  8p 
a  ^P 


n 


TP 

p 


(i) 


(2) 

(3) 


(4) 


(5) 


(6) 


(7) 
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Consider  the  right-hand  side  of  (4): 


£  (PU)  +  F7  (pv) 


3p  9p  ,  /3u  ,  3v\ 

u  53  +  v  +  P(W  +  By) 


From  (1) 


3u  +  3v  (  u  3p  ^  v  3p 


^  By  '  PW 

and  with  this  relation  one  can  write: 


4  (PU)  +  *y  (PV) 


„  9P  ,  v  9P  Pu  9P  Pv  9P 

U^  +  v  T^y 


+  v^py  • 


9 


Hence,  (4)  may  be  written  as 

9  a 

u  q—  +  v  q— 
dx  dy 

_  T  8  3  d 

Now  u  BSE  v  By  "  dS 

where  S  is  the  arc  length  along  the  particle  path.  Therefore,  (4)  becomes 

tj^  n 

e  +  -g-  +  p  "  constant  along  the  particle  path.  (8) 

The  first  law  of  thermodynamics  may  be  written  as 


Tdr)  -  dc  +  pd  i 


(9) 


with  T  =  absolute  thermodynamic  temperature.  In  terms  of  derivatives  along 
the  particle  path  (5)  is 

p  £L?2  o 
p  dS  -  ' 

and  if  one  considers  a  system  in  which  dp  =  0  along  the  particle  path,  then  (9) 
becomes 


From  (8): 


de  =  -pd  —■  along  path  lines. 


dc  +  u  du  +  v  dv  +  i  dp  +  pd  ^  =  0  along  path  lines. 


(10) 


(11) 


Using  (10)  and  (11)  yields 


u  du  +  v  dv  =  -  —  dp  along  path  lines. 


(12) 
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Since  (12)  is  valid  along  path  lines,  one  can  also  write 

du  .  dv  1  dp 
udS  +  vd£=  -pH3* 

i0  0 

Remembering  that  =  u  +  v  (13)  becomes 


2  3u  ,  3u  ,  3v  ,  2  3v  _  1  dp 

u  Bi+UVB^  UVB^+V  By  ~  p  dS  ’ 


and  since  along  path  lines  dp  =  0,  (7)  may  be  written  as 

=  IT"  -  along  path  lines 
P  T)  P 

or,  what  is  the  same  thing,  as 

dp  =  2  dp 

dS  dS  - 


By  using  (15),  (14)  becomes 


2  9u  /3u  3v\  ,  2  3v  a  dp 

u  B3c+uvVBy  BxJ  v  ~  p  dS  * 


Now  (1)  may  be  written  as  follows: 

dp  8u  ,  3p  ,  9v  _  dp  /3u  ,  3v\  n 
uB3i+pB5Z+vByp37~dS  *\Bx  Byj  "  0 


which,  solved  for  yields 


dp  _  /3u  ,  3v\ 

dS  =  ”  P  \Bx  +  Byj  • 

Putting  this  last  expression  in  (16): 


2  fdu  ,  3v\  2  3u 

a  (.55  +  sy)  =  u 


3u  ,  3v\  ,  2  3v 


uv  Uv  +  BxJ  +  v 


/  2  2*  3u  /3u  ,  3v\  ,  ,  2  2v  3v  _  n 

(a  ~  u  ^  Bx  ~  uv\By  Bxj  (  ~  ^y  °  ‘ 

If  one  considers  (6)  along  path  lines  (dp  =  0),  then 

P.=(JL)7  orp  =  pn(-E-ffy, 

Po  'p0  °VV 
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where  pn,  pn  are  the  values  at  some  initial  point,  and  using  (18)  in  (7): 

U  U  ^  ^ 

2-TP/pth  =a2fP-')T~ 

"  ^  W  Po  W  0  W  ’ 


whence 


/a2\T_1  h  2  TP0 

P  =  Pol  —  I  -  where  a, 


0\  2 
-a0 


0  p 


0 


(19) 


Operating  on  (19): 


7 

T_1 


dp  “  P 


7 


da 


0  y  -  1 


7  da 
— —  P 


2 


2  7-1^2 
a0  a 


or 


dp  7  da 


From  (7): 


1  ,  a  , 

—  dp  = - dp  , 

P  7P 


and  using  (20): 


dp 


p  *  7-1 

Putting  this  last  in  (12): 


1  a  2 
da 


2 


1  ,,  2  2,  da 

2  d(u  +  v  >  =  - 


which,  when  integrated  along  path  lines,  yields 

7 


2  2  2a2 

u  +  v  = - j  +  constant,  along  path  lines. 


(20) 


(21) 


For  use  elsewhere,  note  that  if  7  =  3: 

2  2  2 

u  +  v“  +  a  =  constant,  along  path  lines. 


(22) 


2  2  2  2 

Note,  also,  that  (21)  implies  that  a  =  a  (u  +  v  ),  so  that  (17)  is  an  equation 
in  (u,  v)  as  functions  of  (x,  y).  In  order  to  use  the  method  of  characteristics 
here  we  need  another  differential  equation  in  (u,  v)  as  functions  of  (x,  y). 

To  that  end,  consider  the  conservation  equations  for  a  general  hydro- 
dynamic  material  with  t1^  =  -  p 6^  and  f1  =  0.  First  convert  these  equations 
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from  tensor  to  vector  notation  by  noting  that: 

and 

OV,  "  -(P6iju  ).j  ■  -(pu)  •  -V  •  <pu>  ■ 

J  i  o  3  3 

With  these,  the  conservation  equations  become: 
p  +  PV  •  u  =  0, 
pu  +  Vp  = 

pc  +  p  ~2  =  -  V  •  (pu), 

pfj  >  0, 

where  (  )  1  -jjj-  -  -jjj-  +  u  *  V  for  unsteady  motion, 

and  (  )  clg  ~  ’  7  lor  steady  motion. 

For  steady  motion,  (23)  may  be  written  as 
u  •  Vp  +  pV  *  i*  =  0 
or 

7  •  u  =  u  •  Vp  • 

Considering  the  right-hand  side  of  (25), 

V  •  (pu)  =  pV  •  u  +  (Vp)  •  u  , 
and  using  (27), 

V  •  (pu)  =  u  •  Vp  +  jj  •  Vp  =  p«  •  (^Vp  "  "r  Vp) 

or 

7  *  (pu)  =  pu  •  (v  • 

Using  (28)  in  steady  (25): 

Pi!  *  7  (e  +  IT  +  f)  =  0  • 


(23) 

(24) 

(25) 

(26) 


(27) 


(28) 


(29) 
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Now  (9)  may  be  written  as: 

T  ®  =  a!  +  P  ®(;5-)  or  u  -  TVn  =u  -  Ve  -  pVi  . 

From  (24)  for  steady  motion, 

u  =  (u  •  V )  u  =  "  \  V  P> 

r 

and,  in  general, 

V(o+P=V€+ivp  +  PVi. 

so  that. 


+  u  •  (Ve  +  PV^  +  ~  Vp)  • 


Using  (30)  in  this: 

2  2 

U  •  V  (e  +  \  +  g)  =  U  .  V  \  +  u  •  (ye  +  py^-)  -  U  •  (u  .  SJ)u  . 

In  general, 

2  u 

(u  •  V)  -  2  u  •  (u  •  y)J  =  u  •  (u  •  V)u  , 

so  that  (32)  becomes 

“•  V  (5  +  ^  +  £)  =  u  .  (y«  +Pyi)  • 

Using  (30)  in  (33)  and  recalling  (29): 

2  2 

u  •  Ty n  =  u  •  y(e  +  ^  (e  +  \  =  0,  along  path  lines. 

*  j 

Since  u  •  Tyrj  =  T  u  •  \jr\  -  T  n,  (34)  implies: 

2 

e  +  ^  =  constant  and  n  -  constant  along  path  lines. 


(30) 


(31) 


(32) 


(33) 


(34) 


(35) 


This  is  the  same  result  obtained  in  (8);  it  was  repeated  because  the  vector 
equations  were  needed  for  what  follows. 
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Now  the  expression  (u  *.  \J  )u  while  in  vector  notation  is  a  tensor  the  way 
it  stands.  In  order  to  write  it  as  a  vector  expression,  we  revert  to  tensor 
notation  for  a  moment  to  find  that 


(u  •  V  )u  =  u 


i  8(u1s;l) 


ax1 


1  ] 

u  u.  .  eJ 


(36) 


where 

u1  s  contravariant  components  of  u, 
in  a  covariant  components  of  u, 

e"1  =  reciprocal  base  vector. 
Operating  on  the  right-hand  side  of  (36): 


*  •  ■  *  *  *  •  ■  1  O 

1  u.  .  e^  =  u1u.  .e^  -  u\i.  .e^*  +  ill.  .e^  =  Vu 

j;i  ~  i;j~  h.r  r>'~  2 

u1u.  .e^  -  ul.  .e1 

L  i;j~  bj-  J 


1  „  2 

=  9  VU 


_  ktf  i 
6  ..  u  u 
Ji 


St,  k~ 


(see  footnote*) 


l  ki  m  ,] 

e  •  ■  u  e  u  n  i  e J 
jim  i;  k~ 


Now  e1^mu^  ^gj  =  curl  u  =  w  =  vorticity  vector,  and  ej^uV^e*1  - 
where  A  5  cross  product  or  vector  product,  hence 


u  A  w 


6^  is  the  generalized  Kronecker  symbol  which  has  value  1  for  j  -  k,  i  -St, 
and  -1  for  j  -  SL ,  i  =  k. 

...  and  ejlk  (not  to  be  confused  with  e  for  energy  )  are  related  to  the 


jik 

permutation  symbols  e 

r 


kim  P 
e  as  follows 


lif  m  , 

e  ’  kirn  { 


kim’ 

=  1  when  kim  is  an  even 
permutation  of  123 

=  -1  when  kiin  is  an  odd 
permutation  of  123 

=  0  otherwise 


N 


andEjik  =  ejik'S 


.jik  =  e 


jik 


fg 


j 


where  g  =  determinant  of 
the  metric  tensor. 
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Therefore, 


*  *  1  2  1  2 
uV.jg'1  =  7j\/  u  -  uAu  :  gV11  -  H  A  curl  U  • 

(36)  becomes 


1  2 

(u  '  SJ)u-  ■ j  Vu  "  iJ  A  curl  u. 

Fi'om  (30) 

TV»  =  V'  +  pv|  • 

Prom  (31)  and  (38) 


(37) 

(38) 


12  1 
2  V  11  “  i!  A  curl  u  =  “  p  V  P* 

hence 

v(^+f)=Vc  +  PVf  +  jVP 

and 

i\/u2  -  u  A  curl  u  +  V(e  +  “)  =  ?e  +  =  Tyq 

or  2 

TVn  =  V  (e  +  -^-  +  jj)  -  u  A  curl  ~  •  (39) 


From  (39)  it  is  easily  seen  that  if  yp  =  0  everywhere  and  curl  u  =  0, 


then 


everywhere 


or 


c 


+ 


2 

~  +  —  =  constant, 

2  p 


throughout  the  material. 


Therefore',  for  an  isentropic,  irrotational,  steady  process  e  +  u*"7  2  +  p/p 
=  constant  everywhere  and  the  relations  (8)  through  (22)  are  seen  to  hold 
everywhere  instead  of  just  along  path  lines  with  the  additional  requirement 
on  the  (u,  v)*  s  that 


-  curl  u 


3u  3v 
’  Ex 


(40) 


liquations  (17)  and  (40)  are  seen  to  be  the  desired  pair  of  equations. 
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APPENDIX  2 

DERIVATION  OF  THE  CHARACTERISTIC  EQUATIONS'" 

The  set  of  two  differential  equations  to  be  solved  in  this  instance  are; 


(a2  -  u2>  fj  -  uv  ||  -  uv  |H  +  (a2  -  v2>  ||  -  0, 

8u  +  av_8u+08v=0. 

ox  ox  o  y  dy 


(i) 


In  general,  a  set  of  first-order  partial  differential  equations  in  two  independent 
variables  may  be  written  as 

(2) 


L[  U]  =  C 


where 


L[  U]  -  AUx  +  BUy  , 

A,  B  =  matrices, 

C  =  vector, 

U  =  vector  representing  the  dependent  variables, 

9U 

It  is  readily  seen  that  to  write  (1)  in  the  form  of  (2),  one  defines 
U  =  (u,  v]  , 


A  = 


-uv 
1  . 


B  = 


-uv 

-1 


0 


C  =  [  0  0]  . 

In  this  case  (2)  is  of  the  form 

AU  +  BU  =  0, 

—  x  ~  y 

which  may  be  written  as 

U  +  A"1  BU  =  0 

~x  ~  y 


(3) 

(4) 


The  treatment  here  follows  closely  that  in  Ch. 


5  of  ref.  7. 
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))rovided  |  A  |  0,  where 


I  a|  =  determinant  of  A, 


A  ^  -  inverse  of  A. 

Suppose  it  is  given  that  U  has  some  known  value  on  a  curve 
C[  4>  (x,  y)  =  0]  along  which  the  inner  (or  tangential)  derivative  (see  ref.  7,  pp. 
132-133)  of  U  is  also  known.  In  that  case,  the  outward  normal  to  C  is  given 
by 


V<j>  = 


3<j>  .  . 

FyJ 


(5) 


and  the  unit  normal  vector  is  given  by 


The  derivative  of  a  function  F  along  a  curve  with  S  the  arc  length  is 
given  by 


+ 


9F 

y  ~ 


3v  1 


0 


s 


where  S  is  the  unit  tangent  vector.  If  S 
known  that 


(7) 

[S i ,  S2]  and  n  =  [n^,  n2j,  it  is  well 


■n2  and  n^  =  S2; 


therefore  if  we  let  -  <)>.,  etc.  then 


S  = 


so  that  (7)  becomes 


F  <j>  +  F  <j> 
x  y  y  x 

J<t>  2  +  4>  2 

v  x  y 


(8) 


Since  we  are  assuming  the  inner  derivative  of  U  is  known  on  C,  (8) 

2  2 

used  to  express  the  inner  derivative  of  U  assuming  <j>x  +  f  0: 

-  U  <t>  +  U  <j>  =  known  quantities, 

y  ~y  x 


can  be 

(9) 
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and  since  cj>(x,  y)  is  known,  4>  and  4>  are  known,  so  define  . 

x  y 

*x 

T  r 

•  y 


(10) 


and  using  (10)  in  (9).  and  solving  for  Ux, 

U  -•  -t  U  +  (known  quantities),  and  substituting  (11)  in  (5),  (11) 

~  x  ~  y 


■  rU  +  A  1BU  '+  (known  quantities)  =  0. 

y  y 


(12) 


In  order  to  determine  U  from  this  expression,  the  determinant  Q  of  the  co¬ 
efficients  of  Uy  must  satisfy 

Q  =  |a-1B  -  tI  f  0,  where  I  =  identity  matrix. 

If,  on  the  other  hand,  Q  =  0,  the  curve  C  is  called  characteristic  and  the 
characteristic  condition  is 


Q  =  |a_1B  -  t I 


=  0. 


(13) 


The  system  is  classified  as  hyperbolic  if  (13)  has  two  distinct  real  roots.  In 
that  case  a  solution  by  the  use  of  characteristics  is  possible.  Using  (3): 


| A  |  =  a2  -  u  , 
1 


A 


-1 


uv 


2  2  2  2 
a  -u  a  -u 


A_1B  = 


0 

-2uv 


2  2 
a  -v 


2  2  2  2 
a  -u  a  -u 


.,2,2 

if  a  f  u  , 


(14) 


Putting  (14)  in  (1 3): 


Q  = 


-  o  2  2 

,  -2uv  a  -v 

~2  2  ~2  2 
a  -u  a  -u 


i 

'  l 

0 

-  T 

.0 

l . 

=  0 


-1  0 

The  characteristic  condition  (15)  is  solved  as  follows: 


(15) 
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From  (16)  it  is  clear  that  there  are  two  distinct  real  roots  if 

2  2  2  n 
u  +  v  -  a  >  0  . 

If |a ”^B  -  Tl|=  0,  then  there  exists  a  vector  SL  such  that 
iA -  frl  =  0  or  IA  =  frl  . 


Letting  SLa  correspond  to  t  (o-  =  1,  2),  (18)  may  be  written  as 

f°(A_1B)..  =  ftaT  6  ••  =  £ar 
1  i]  1  O'  ij  J  O' 

from  which  one  must  find  Sta .  Expanding  (19)  and  solving: 

‘l  [CA"lB)n  -TJ+,2(A"1b)21  =  0' 
<(A'lB)l2+i2[(A‘lB)22-TJ 


whence 


(a-1b) 
a-'b).. 


(a-'b). 


21 _ £<*  Aa  =  v _ hi-  T_g 

11  "T*  CA  BA2 


Using  (14)  and  (16)  in  the  second  of  (20): 


1,2  _  _  uvt  a/u2  +  v2  -  a2  1,  2 

^1  2  2  2  2  ^  2 
a  -  v  a  -  v 

2  2 
a  -  u 

Now,  operate  on  (5)  with  and  make  use  of  (18): 

jP.^U  1  +  j?^(A -1b).  .U**  =  f^U1  +  Jta  t  6..U  j 
ix  i  x  7  rj  y  ix  i  o  l]  y 


I^U1  +  Jt ar  U1  -  0  . 
ix  l  a  y 


(16) 


(17) 


(18) 


(19) 


(20) 


(21) 
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This  last  may  be  written  as 
ri 


(22) 


the  derivative  along  the  H01  direction. 

Using  (21)  in  (22): 

uv  ~  a  2+  v2g~  —  D^1  +  D^2  =  0  along  1, 
a  -  v 


(23) 


and 


,  /  2 
uv  +  avu 


2  2 
a  -  v 


D 


2UX  +  D2U2  =  0  along  2. 


In  order  to  determine  the  curve  C,  think  of  §(x,  y)  -  0  solved  for 
y  =  y(x),  and  since 

d6  =  ||  dx  +  ||  dy  -  0, 


3y 


ais  -q;  =  T  £see  <10)] 


(24) 


Using  (16),  (24),  and  (23)  one  can  then  write  the  differential  equations  for  the 
characteristic  curves  and  the  compatibility  equations  which  must  be  satisfied 
along  them  as: 

2  2 
v  u 


r  dy  _  _  uv  +  avu  +  v  -  a 

C2:  dx  T2 


2  2 
u  -  a 


along  which 


uv 


+  as/u2  +  v2  -  a2  du  dv 


2  2 
v  -  a 


dx2  dx2  ’ 


(25) 


and 


/  2  2  2 

dy  _  ,  _  uv  -  avu  +  v  -  a 

Cl:  d5E  T1  2  2 

u  -  a 


uv 


/  2  .  2 
avu  +  v 


2  2 
v  -  a 


du 

dx. 


dv 

dx^  ' 


along  which 
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APPENDIX  3 

DERIVATION  OF  THE  CHARACTERISTIC  EQUATIONS 
IN  THE  HODOGRAPH  PLANE  ' 


The  set  of  partial  differential  equations 

3v  9u  ,  /_  2  __2X  9v 


j  2  2v  C/U  ^  V  V/  l  \  I  Lj  \ 

(a  -  u  )  -  uv  ^  -  uv  ^  +  (a  -  v  ) 


>y 


3y 


=  0. 


(1) 


9u  9v  9U  Q  9v  =  Q 

°'97Z  F^*^  0'9y  °» 


is  a  set  of  equations  to  be  solved  for  u  =  u(x,  y)  and  v  =  v(x,  y).  Since  there 
are  two  dependent  and  two  independent  variables  it  is  possible  to  transform 
(1)  into  a  set  which  is  to  be  solved  for  x  =  x(u,  v)  and  y  =  y(u,  v).  To  that  end, 
we  write: 


,  9u  ,  ,  9u  , 

du  ,  dx  +  dy  , 

,  9v  ,  9v  , 

dv  -  -js-tt  dx  +  -g— ;  dy  , 


(2) 


Fx  tfy 

and  consider  (2)  as  a  set  of  equations  to  be  solved  for  dx  and  dy.  Hence 


dx  = 


dy  = 


du 

9u 

dv 

9  v 

^y 

j 

9u 

Fx 

du 

9v 

FV 

dv 

1  9v  ,  1  9u  , 

J  Fy  dU  ~  J  Fy  dV' 


(3) 


1  9v  ,  ,  1  9u  , 

j^du  +  733dv- 


where 


3u  9v  9u  9v  ,  n 

Fx  9y  9y  9x  ' 

By  considering  x  =  x(u,  v)  and  y  =  y(u,  v),  one  can  also  write 

,  9x  .  9x  . 
dx  =  g-jj  du  +  dv. 


(4) 


_  9y 


9y 


dy  =  du  +  -75^-  dv. 
J  du  dy 


This  treatment  also  follows  ref.  7. 
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Equating  like  terms  in  (3)  and  (4)  and  assuming  J  f  0: 

9x  _  9v  T  9x  _  9u 
J  Fu  ~  Fy’  J  Fv  ~  ~  Fy’ 

3y  9v  3y  _  3u 
J  Fu  ^  ~  Fx’  J  Fv  ~  Fx  ' 


Putting  (5)  in  (1)  and  rearranging  terms: 

,2  2.  9x  ,  3y  3x  ,  ,  2  2^  3y  _  n 

(a  -  v  )  ^  +  uv  ^  +  uv  +  (a  -  u  )  -  0, 


Equations  (6)  can  now  be  treated  just  as  (1)  was  treated  in  appendix  2. 


z  =  [x,  y]  , 

°i  -  FT’ 

r  2  2 

a  -  v 

uv 

A  = 

0 

-1 

r  2 

uv  a 

2 

-  u 

B  = 

and  (6)  becomes 


AZ  +  BZ  =  0  . 
~  u  ~  v 


Equation  (8)  can  be  solved  for  Z  which  yields 


Z  +  A  BZ  =  0 
~  u  ~  v 


where 


A-1B 


o  2  2 

2uv  a  -  u 

2  2  2  2 

a  -  v  a  -  v 


if  |  Aj  =  -  a?  f  0.  & 


L  -1  0 

The  characteristic  condition  is: 


2  2 
a  -  u 


i  2  2 

Q  =  A  !B  -  tI  =  a  V 


2  2uv 


"  “  y  rr  4-  “ 

~2  2  T  T 

a  -  v  a 
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which  has  roots 


1,2 


JI 


2  2 
uv  ±  ax/u"  +  v  -  a 


2  2 
a  -  v 


This  system  is  therefore  hyperbolic  if 


2  2  2  .  n 
u  +  v  -  a  >  0. 


The  null  vectors  £a  are  then  found  to  be: 


f1'2 

X1 


J  2  ,  2  2  1  0 

uv  ±  a^u  +  v  -  a  „1,2 

2  2  2 
a  -  u 


Using  Sa  in  (9),  the  compatibility  equations  are; 


'  ^  J  2  T  2  2\ 

uv  +  ax/u  +  v  -_a_jDix=  Dly, 


2  2 
u  -  a 


v/  2  ,  2  2 
uv  -  avu  +  v  -  a 

2  2 
u  -  a 


^2X  ^2^" 


(12) 


(13) 


(14) 


(15) 


If  the  characteristic  curves  in  the  u,  v  plane  are  designated  by  T,  then: 

/uv  +  a/u2" 


r  •  =  t 

r  du  i 


,  2  2' 
+  v  -  a 


~~2  T 
v  -  a 


along  which 


uv  +  avu 


u  -  a 


dx 

du. 


dy 

du. 


and 


J? 


avu  +  v 


along  which 


avu  +  v 


2  2 
u  -  a 


dx 

du2 


dy 

du2  ' 


(16) 


A  comparison  of  equations  (16)  above  and  (25)  in  appendix  2  reveals  the  rela¬ 
tion  between  the  characteristic  curves  in  the  x,  y  plane  and  the  u,  v  plane.  F 
curves  in  the  u,  v  plane  are  compatibility  conditions  in  the  x,  y  plane  and  C 
curves  in  the  x,  y  plane  are  compatibility  conditions  in  the  u,  v  plane. 


-68- 


APPENDIX  4 

CONDITIONS  AT  THE  DETONATION  FRONT 

Consider  a  detonation  front  traveling  in  the  -x  direction  through  a 
rectangular  slab  of  material  whose  z  and  x  dimensions  are  infinite.  Det  D  be 
the  speed  with  which  the  front  moves.  To  consider  this  as  a  steady  state 
problem,  the  detonation  front  is  brought  to  rest  and  assumed  to  be  located  at 
x  =  0.  The  following  diagrams  may  prove  helpful  in  visualizing  the  process: 


y 


Fig.  20 

The  transformation  -u'  +  D  is  made  to  bring  the  front  to  rest  and 
(x'  -  Dtf)  is  made  to  put  the  front  at  the  x  =  0  position.  Since  the  slab  is 
infinite  in  the  x  direction  this  last  does  not  really  play  a  part.  The  diagram 
now  becomes: 


y 


Fig.  21 
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At  the  front  now  the  particle  velocity  of  burnt  material  is  u  =  D  -  u'  in  the  x 
direction,  v  =  0  in  the  y  direction,  where  u  is  positive  if  it  is  in  the  positive 
x  direction  and  v  is  positive  if  it  is  in  the  positive  y  direction. 

The  problem  considered  in  this  paper  deals  with  the  material  behind  the 
front  (burnt  material)  and  this  latter  coordinate  system  is  the  one  used.  In 
that  case  the  velocity  u  is  the  velocity  in  the  x'  direction  relative  to  the  detona¬ 
tion  front  so  that  at  x  =  0,  u  =  a  .  where  a  .  is  the  sound  speed  at  the  detona- 

C_J  C-J 

tion  front  as  determined  by  the  Chapman- Jouguet  hypothesis  (see  ref.  2,  p. 


212). 


In  order  to  determine  a  .  and,  therefore,  u  at  x  =  0  consider  Fig.  20 

^  ~  J 

above  in  conjunction  with  the  Rankine-Hugoniot  equations,  with  the  detona¬ 
tion  speed  D  in  place  of  the  shock  speed  (see  ref.  1,  p.  3),  that  hold  at  the 
front: 


(u'  -  Uq)p 

P  "  P0 


(1) 


p  -  p0  =  p0d(u'  -  uo)  • 


(2) 


For  an  ideal  gas  the  sound  speed  at  the  front  is  given  by 

2  y  p 
a  .  =  . 

c-3  P 


(3) 


The  statement  u  =  a^  ^  is  equivalent  to  the  Chapman- Jouguet  hypothesis  that 
holds  at  the  front  and  may  be  written  as 


a  .  +  u'  =  D. 
c-3 


(4) 


By  considering  the  unburnt  material  to  be  at  rest  with  p^  «  p,  equations 
(1)  and  (2)  become: 


y 

(5) 

9  ~  P0 

=  P0Du»  . 

(6) 

From  (5)  solve  for 

^0  _  D  -  u» 

P 


D 


(7) 
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From  (4)  solve  for 


a  .  =  D  -  u'  . 

c-] 

Using  (3)  in  this  last  expression: 

(D  -  u'  )2  =  —  • 

P 

Putting  (6)  in  this  expression: 

o  TPnDu' 

cd  -  u«  r  =  -4 — 


which  may  be  solved  for 


P_0  _  (D  -  u »r 

p  y  Du' 

Combining  (7)  and  (8)  yields 

t 

D  -  u'  _  (D  -  u')‘ 


(8) 


D 


which  becomes 


y  Du' 


D  -  u'  _  1 

yu'  ' 

If  now  this  expression  is  solved  for, 

D 


u' 


+  1 


(9) 


and  u'  is  substituted  from  (9)  into  (8),  a  .  may  be  written  in  terms  of  D  as 

c_  J 


-D. 


and 


c-j  y  +1 

Hence,  the  boundary  conditions  at  the  detonation  front  are 

u  -  a  .  =  X  i  D  at  x  =  0 
c-3  7  +  1 

v  =  0  at  x  =  0. 


(10) 


(ID 


If  7  =  3,  (11)  is  seen  to  be 

u  =  a  .  4  D  at  x  =  0, 
c-J  4 

v  =  0  at  x  =  0, 


(12) 
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APPENDIX  5 

GENERAL  RELATIONS  BETWEEN  THE  CHARACTERISTIC  CURVES 
IN  THE  x,  y  PLANE  AND  THOSE  IN  THE  u,  v  PLANE  ' 

Introducing  parameters  a,  j3  along  characteristic  curves  C+  and  C 
respectively,  the  characteristic  equations  (25),  appendix  2,  and  (16),  appendix 
3,  may  be  written  as: 


C+:  y 


a 


J2  ,  2  _ 2 

Vu 


uv  +  avu  +  v 

2  2 
u  -  a 


x  , 
a 


r+:  v 


a 


uv  +  a 


[2 

vu 


,  2  2 
+  v  -  a 


2  2 
v  -  a 


u  , 
a 


c_:  yp  - 


'  .f  2  ,  2  2 

uv  -  avu  +  v  -  a 

2  2 
V  u  -  a 


(i) 


r":  V 


uv  -  a 


J2 

vu 


.  2  2 
+  v  -  a 


0 


2  2 
v  -  a 


v 


where  A 


dA 


,  etc. 


a  do- 
Using  (1),  one  can  form 


u  xn  +  v  yQ  =  u  Xn 
a  p  cr  p  a  p 


,  /  2  .  2  2  V  2  ,  2  2 

uv  +  avu  +  v  -  a  uv  -  avu  +  v  -  a 


2  2 
v  -  a 


2  2 
u  -  a 


=  u  x„ 
a  p 


22  2/2.2  „2v 

uv  -  a(u  +  v  -  a) 


22  22  2  2.4 

uv  -  au  -  av  +a 


=  0, 


or 


u 


a 


a 


x 


£ 


which  may  be  written 


du 

dv 


a 


(D 


0 


(2) 


From  (2)  one  concludes  that  if  T  and  C  are  plotted  in  the  same  co¬ 
ordinate  system,  then  F  is  perpendicular  to  C  . 


The  discussion  here  follows  closely  that  in  ref.  2,  p.  259  ff. 
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Using  (1),  one  can  also  form 


u0x  +  v0y  =  u0x 

(3  a  (3Ja  (3  a 


uv  -  a 


/~2  2  2  ,  /  2  ,  2  2 
Yu  +  v  -  a  uv  +  av  u  +  v  -  a 


2  2 
v  -  a 


2  2 
u  -  a 


=  0, 

so  that  one  can  write 

§2)  = _ I _  (3) 

4  (§)i; 

From  (3),  one  concludes  that  if  F"  and  C+  are  plotted  in  the  same  co¬ 
ordinate  system,  then  T"  is  perpendicular  to  C  . 

Consider  a  point  x,  y  in  the  x,  y  plane.  The  flow  direction  0  is  given  by 


tan  0  =  —  .  (4) 

Let  the  angle  between  the  flow  direction  at  x,  y  and  the  tangent  to  the  C  curve 
at  x,  y  be  A,  and  the  angle  between  the  horizontal  and  the  tangent  to  the  C 
curve  be  <j>+,  then 

c|>+  =  0  +  A  (see  Fig.  22),  (5) 


whence 


tan  <j>+  =  tan  (0  +  A)  = 


tan  0  +  tan  A 
1  -  tan  0  tan  A 


(6) 


Let  the  angle  between  the  flow  direction  at  x,  y  and  the  tangent  to  the  C  curve 
at  x,  y  be  B,  and  the  angle  between  the  horizontal  and  the  tangent  to  the  C 
curve  be  ,  then 


<|>~  =  0  -  B  (see  Fig.  22), 


(7) 


whence 


tan  0  -  tan  B 
1  +  tan  0  tan  B 


tan  <j>  =  tan(0  -  B)  = 


(8) 
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Fig.  22 


From  equation  (1)  above  it  is  clear  that 


y 


tan  <t>  = 


+  a  uv 


/  2  2  2 
’  +  avu  +  v  -  a“ 


a 


2  ,2 
u  -  a 


and 


tan  4>  = 


yQ  /2T2  2 

J  p  uv  -  avu  +  v  -  a 


x. 


2  2 
u  -  a 


From  a  comparison  of  (9)  with  equations  (15)  and  (16)  in  appendix  2: 


9  2 

4  2  . +  ,  2uv  .  .  +  ,  a“  -  v  _  n 

tan  4>  +  — o  tan  <p  +  — ^ ^  -  0, 


2  2 
a  -  u 


2  2 
a  -  u 


and 


tan^  4> 


o  2  2 

2uv  ,  ,  -  ,  a  -  v  _  n 

tan  0  -r  — , - ?r  -  0. 


2  9 

a  -  u" 


2  2 
a  -  u 


Using  (6)  in  (10. a), 


9 


(9) 


(10) 


/  2  2,  /tan  0  -  tan  A  V  tan  9  +  tan  A  2  T.^  _  n 

(a  -  u  )  ( - - —  5  —  \  )  -  -uv  - - t  +  a  -  v  -  0. 


^1  -  tan  9  tan  Ay 


1  -  tan  0  tan  A 


Clearing  fractions  and  squaring 


999  99999  2  9 

a“  taiCO  +  2a-  tan  0  tan  A  -  a"  tan  A  -  u“  tan-@  -  2u“  tan  9  tan  A  -  u  tan- A 


2  2  2 
+  2uv  tan  9  +  2Uv  tan  A  -  2uv  tan  9  tan  A  -  2uv  tan  9  tan-A  -  a“ 


-  2a^  tan  9  tan  A  +  tan^O  tan^A  -  v^  +  2v“  tan  9  tan  A  -  tan^O  tan^A  =  0 


-74- 


Using  (4)  in  the  above  after  eliminating  the  terms  whose  sum  is  zero; 
2 

a2  +  a2  tan2A  -  v2  -  2uv  tan  A  -  u2  tan2A  +  2v2  +  2uv  tan  A 


u 


3  2  3  4 

-  2  tan  A  -  2v^  tan^A  +  +  a^  tan^A  -  v  +2  —*  tan  A - ^  tan  A  =  0, 

u  z  u  * 


u 


u 


from  which 


2 

^a2  +  a2  (1  +  tan2A)  =  tan2  A 


2  n  2  v 
u  +  2v  +  — 

c 

u 


Solving  this  last  for  a  : 
2 


tan2A 


1  +  tan  A 


u2  +  2v2  +  y4/u2  =  ( tan2 A  \  u4  +  2u2v2  +  v4 


1  ,  2/2 
1  +  v  /u 


V 


2a 
sec  A 


2,2 

u  +  v 


2  2  2 
and  letting  q  =  u  +  v  , 

2  2.2.  /i  i  \ 

a  =  q  sin  A  .  U1J 

Using  (4)  and  (8)  in  (10. b)  by  a  similar  calculation: 

a2  =  q2  sin2B  .  (12) 

Now  at  the  point  (x,  y),  a  and  q  have  some  fixed  value  so  that  except 

2  2  v 

where  a  =  u,u  =  0orl±  —  tan  A  =  0, 

A  =  B.  (13) 

Hence,  with  the  exceptions  noted,  the  flow  direction  bisects  the  angle 
between  C+  and  C”  and  since  C+  is  perpendicular  to  T"  and  C”  is  perpendicular 
to  r+,  the  flow  direction  bisects  the  angle  between  F  and  F  . 

Let  the  angle  between  the  flow  direction  and  the  F+  curve  at  the  point 
(u,  v)  in  the  u,  v  plane  corresponding  to  the  point  (x,  y)  in  the  (x,  y)  plane  be 
A*',  and  the  angle  between  the  horizontal  and  the  tangent  to  the  T  curve  be 
i[j+ ,  then 

ijt  =  e  +  A' 

or 

+  _  tan  0  +  tan  A1 
1  -  tan  6  tan  A^ 


tan  [p 


(14) 
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From  equation  (1) 


tan  tp  = 


v  /  ,  /  2  ,  2  2 

^  /  uv  +  avu  +  v  -  a 


u 


a 


2  2 
v  -  a 


and  comparing  (15)  with  equations  (11)  and  (12),  appendix  3, 


(15) 


,  2,+  2uv  ,+  a  -  u  _  n 

tan  ip - ry - 2  tan  ^  +  ~2 - 2  _  U* 

&  -  v  a  -  v 


(16) 


Using  (14)  and  (4)  in  (16)  and  solving  for  a  : 


1  +  (1  +  tan^A1)  =  u^  +  2v^  +  ^ 

uv  u 


or 


2  2  2  A 
a  =  q  cos  A1 


(17) 


(which  is  another  proof  that  F  curves  are  normal  to  C  curves).  Comparing 
(17)  and  (1 1), 


A-  =  |  -  A, 


(18) 


and  one  concludes  that  the  component  of  flow  normal  to  a  C  curve  is  equal  to 
the  sound  speed  and  the  component  of  flow  tangent  to  a  F  curve  is  equal  to  the 
sound  speed. 

From  (2),  (3),  (9),  (14),  and  (15)  it  is  clear  that  (1)  can  be  written  as 


.+ 


+ 


r  :  v  =  tan  ip  u  , 
a  a 


C+:  y 


x 


a  ,  ,  -  a 

tan  ip 


r”:  Vp  =  tan  ip  u ^  , 


C":  y 


0  tan  0+  P  ' 


(19) 
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APPENDIX  6 

ANALYTIC  SOLUTION  FOR  CHARACTERISTICS 
IN  THE  HODOGRAPH  PLANE" 

To  obtain  the  solution  for  T  curves  in  the  u,  v  plane,  start  with  equa¬ 


tions  (15)  and  (17)  from  appendix  5: 


Along  r+:  ^  =  tan  ip+. 

(1) 

At  (u,  v):  cos^  A' . 

(2) 

Again  referring  to  appendix  5: 

ip+  =  9  +  A'  =  0  -  A  +  J  . 

(3) 

From  (1): 

du  sin  ip  -  dv  cos  * p  =0. 

(4) 

From  (2)  since  a  >  0,  q  >  0  and  |A'|  <  ^  : 

a  =  q  cos  A'  . 

(5) 

From  (3): 

A'  =  ip+  -  0, 

whence 

cos  A'  =  cos  (i p  -  9)  =  cos  ip  cos  9  +  sin;// 

sin  9.  (6) 

V 

Since  tan  0  =  —  , 

sin  9  =  cos  9  =  where  q  -  vu  +  v  . 

(7) 

Using  (6)  and  (7)  in  (5): 

.  +  4* 

a  =  v  sin  ip  +  u  cos  ip  . 

(8) 

In  appendix  5  it  is  shown  that  a  (s  sound  speed)  is 

the  velocity  along  T 

curves,  and  to  find  the  velocity  g  normal  to  the  T  curves  note  that  the  flow 

velocity  q  is  given  by 

q  =  u  i  +  v  j  , 

(9) 

where  i  and  j  are  unit  vectors  along  u  and  v  respectively,  and  that  the  unit 
vector  along  T+  is 

This  discussion  follows  ref.  2,  p.  264  ff. 
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n  =  i  cos  ip+  +  j  sin  ip+  . 


(10) 


Since  a  is  the  component  of  q  along  r  , 

+  .  + 

a  =  q  .  n  =  u  cos  ip  +  v  sin  ip 


which  checks  with  (8).  To  find  the  normal  component  of  q,  one  takes  the 
inner  product  of  q  with  v ,  the  unit  normal  to  F+.  To  find  v ,  note  that 


dn 


~  dtp 


+ 


-  i  sin  ip+  +  j  cos  i p+ , 


(ID 


whence 


g  =  q  •  v  -  v  cos  i//+  -  u  sin  ip+. 


(12) 


Next,  solve  (8)  and  (12)  for  u  and  v  in  terms  of  a  and  g: 

+  .  + 
a  =  u  cos  ip  +  v  sin  ip  , 

.  4*  4~ 

g  =  -u  sin  \Jj  +  v  cos  ip  . 


The  determinant  of  coefficients,  D,  is 


D  = 


cos  sin  \p  + 

z  +  + 

-sin  ip  cos  ip 


-  1 


whence 


u  = 


a  sin  ip 
g  cos  ip 


4"  .  + 

a  cos  ip  -  g  sin  ip  , 


(13) 


v  = 


COS  Ip 


.+ 


-sin  ip  g 


g  cos  ip+  4-  a  sin  ip+ . 


4- 


(14) 


The  solution,  then,  is  obtained  by  considering  ip  as  the  variable  along 
T+  and  determining  a  and  g  as  functions  of  i p+. 

Using  (13)  and  (14): 


=  _  a  sin  \p+  +  a  +  cos  ip+  -  g  cos  \p+  -  g  +  sin  \p+ , 
dip  \p  ip 

dv  .4-  4-  4-  .  4- 

— ^  =  -  g  sin  ip  +  g  +  cos  ip  +  a  cos  ip  +  a  +  sin  ip  . 

dip  ip  ip 


(15) 


Dividing  (4)  by  dip 


+ 
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du  •  +  dv  ^  +  _  n 

sm  i p - +  cos  p  -  U, 


+ 


di Jj  dip 

which  upon  substitution  from  (15)  yields 

.2+  .+  +  +  +  .  2  .  + 

-  a  sin  ip  +  a  sm  ip  cos  ip  -  g  sm  \p  cos  ip  -  g  +  sm  ip 

ip  p 

+  +  2+  2 ,  +  .  +  ,  +  _  n 
+  g  sin  ip  cos  ip  -  g  cos  ip  -  a  cos  ip  -  a  sm  ip  cos  ip  -  (J 


P 


P 


or 


whence 


2+  2+  2+  .2  + 
a  (sin  ip  +  cos  ip  )  -  g  ,  (cos  ip  +  sin  p  )  -  0, 


P 


P 


+ 


a. 


(16) 


Using  equation  (21)  from  appendix  1  and  equation  (11)  from  appendix  4: 
2  a2  +  u2  +  v2  =  2  «24.„2  =  _2_/a2  W/.2  Wn=l_ti*2 


7-1 


7-1 


2.  2  2  /2\,/2\  f.  y  +  12 

h  7  -  1  V  c-j'  \  c-y  7"1  c-j 


or 


2a  7  +  12  2 

- r  -  - t  3-  .  ~  Q 

7  -  1  7  -  1  c-j 

2  -Y  -  1 

and  if  we  let  p  =  +■  +  p  then 


2  a2  a2  2 

2a2  c-j  2  c-j  2  2a 

—  -  q  or  -V =  q  +7^~i  • 


7  -  1  2 
1  h 


(17) 


Since  a  and  g  are  orthogonal  components  of  q, 

2  2  ,  2 
q  =  a  +  g  , 

so  that  using  (18)  in  (17)  yields 

2  2  2  2 
a.  n2  „a.  o/o,  ,  \  ;  -  a 

c-j  2a  2  _  c-j  ^2(2+y-l\_  c-j 

1~Tl'a  T~  a  V  7  -  1  J  2 

p  '  ju  x  /  iu 


or 


2  2 
ac-j  "  a  2 

- ^ -  =  g  . 

M 


=  g 


(18) 


(19) 
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Differentiating  (19): 
2a  a 


0 


+ 


2g  g  +• 

0 


(20) 


Putting  (16)  in  (20): 


or 


-  a  a  +  =  -p  ga 


a  +  -  H  g 
0 


(21) 


Then  the  problem  is  reduced  to  solving  (16)  and  (21).  In  addition,  note 
the  boundary  condition  from  (19): 


g  =  0  when  a  =  a 


c-3 


(22) 


Operating  on  (16), 


d_J_  +  =  o 

-K2  . .+  ’ 


d  (i/y  )  dijj 
and  using  (21), 


d^g 

d</>2 


whence 


+  M  g  =  0, 


g  =  E  sin  j u(//+  -  (//^,),  where  E  and  ipZ  are  constants 


+ 


(23) 


Using  (23)  in  (21), 


da 


pZE  sin  p(^+  -  ipZ), 


d0' 


+ 


whence 

a  =  -p  E  cos  p(i//+  -  <//l)  +  F,  where  F  is  a  constant. 
Using  the  condition  (22)  in  (23) 


(24) 


-80- 


0  =  g  =  E  sin  fa*  -  i p~l) ,  where  ip*  is  the  value  of  /  at  a  -  ac_^ . 


This  condition  is  satisfied  if 

+  _  ,+ 

i/'q  -  0*  » 


and  using  this  in  (24),  assuming  F  =  0, 

a 

c-J 


a  .  =  -uE  or  E:- 
c-J  M  M 


whence 


ar_.  cos  ju(;//+  -  i/vi)  and  g  =  -  sin  /Lt(^+  -  (//l) 


c-J 

Physically  a  >  0,  so  that 


(25) 


<//+  _  0* 


< 


7T 


2ju  ‘ 


It  is  shown  in  appendix  5  that  0  bisects  the  angle  between  the  T  and  the 
r-  curves  so  that  we  can  write  the  equations  for  T  that  correspond  to  (1), 

(2),  and  (3)  as: 


Along  r  : 


dv 

du 


tan  ip  . 


At  (u,  v):  a2  =  q2  cos2  A', 


From  (26): 


From  (27), 


ip~  =  e  -  A'  . 


du  sin  ip~  -  dv  cos  ip  =  0. 


a  =  q  cos  A'  . 

From  (28), 

A '  =  0  -  v>“, 
and  from  (31), 

cos  A'  =  cos(0  -  i p~)  =  cos  0  cos  0"  +  sin  0  sin  ip  . 


(26) 

(27) 

(28) 

(29) 

(30) 

(31) 

(32) 


Now  by  inspection  it  is  seen  that  (2  9),  (30),  and  (32)  are  the  same  as  (4),  (5), 
and  (6)  with  ip+  replaced  by  ip~ .  Hence,  the  equations  for  the  T~  curves  are 


a  =  a  cos  /l t(0~  -  ipZJ  and  g 
c  3 


c-J 

iU 


sin  fu 


-  t*) 


(33) 
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with 


From  (12)  and  (3),  along  F+: 


+  + 

g  =  v  cos  i //  -  u  sin  \jj  -  v  cos  (0  +  A' )  -  u  sin  (0  +  A' ) 

=  v  cos  0  cos  A'  -  v  sin  0  sin  A'  -  u  sin  0  cos  A'  -  u  cos  0  sin  A'  . 

Using  (7): 

2  2  2 
g  =  —  cos  A'  -  —  sin  A'  -  —  cos  A'  -  —  sin  A'  = - sin  A'  . 

e  q  q  q  q  q 

But  A'  =  ^  -  A  (see  equation  (18),  appendix  5),  whence 

g  =  -q  sin^-g-  -  A)  =  -q  cos  A,  along  T+. 

From  the  dual  of  (12)  for  \jj  and  (28),  along  T  : 

g  =  v  cos  ;//"  -  v  sin  i jj  =  v  cos (0  -  A')  -  u  sin(0  -  A' ) 

=  v  cos  0  cos  A'  +  v  sin  0  sin  A'  -  u  sin  0  cos  A'  +  u  cos  0.  sin  A'  . 

Using  (7): 

2  r  2 
g  =  —  cos  A'  +  —  sin  A'  -  —  cos  A'  +  —  sin  A' 

e  q  q  q  q 

but  A'  =  7r/2  -  A,  whence 

g  =  q  sin  A'  =  q  sin(^  -  a)  =  q  cos  A,  along  T  . 


(34) 


(35) 


From  (34)  and  (35),  since  j  A|  <_  7r/  2,  one  concludes 

g  <  0  for  T+,  and  g  >  0  for  T  curves.  (36) 

Using  (36)  in  (33)  and  (25): 

>0^  an<^  0  .  (37) 

Let  ip  be  the  variable  along  either  F  curve,  then  using  (13)  the  solution  for  T 
curves  is  given  by: 
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u 


c-j 

V 

L 

c-J 


=  cos  n(ip  -  x pj  cos  xp  +  -  sin  q(<//  -  xp  J  sin  i p, 
*v‘  M 


cos  jU  (xp  -  x pj  sin  xp  -  i  sin  ju((//  -  cos  xp. 


(38) 


In  order  to  see  what  the  equations  look  like,  proceed  as  follows: 
— —  =  cos  q(i//  -  i//;J  cos  ip  +  —  sin  n(ip  -  xp  J  sin  xp 

c-J 


J_  ,  1 
2q 


N 

+  cos  ip  cos  n(\p  -  xp  o„)  +  +  ?r)  sin  xp  sin  n(xp  -  xp 

-  (-jj  -  cos  (//  cos  ju(^  -  ip  J  +  -  -J>)  sin  xp  sin  n(xp  -  1//*) 


--2ju  2/ 

_  1/1  ,  ^ 
-  2  W  +  V 


or 

u 


C-J 
Now  let 

then 


U + $ 


cos 


5  =  s(?  -  ac-j  • 


cos 


xp  -  fJ.(xp  -  Xp 


1/1  1 

\  _r 

.  -  2  V  -  1 

)  cos 

xp  +  n{\p  -  xp .) 


((//  -  xp  J(1  -  m)  +  -  ■|('J  -  1)  cos  (1//  -  xp  J(1  +  m)  +  xp,,. 


*6+0 


c-j 


4  f  —  -  l)  a  .  +  a  .  =  q  +  a 
2  Vi  /  c-j  c-j  M  c-j 


(39) 

(40) 

(41) 


Using  (40)  and  (41)  in  (39): 

|~  ( 1  -  ju  )(xp  -  xp  +  xp  ;.J  -  q  c 


u  =  (a  .  +  q)  cos 
c-j 


cos 


(1  +  q)(i p  -  Xp  J  +  xp.. 


From  (40): 


a 


2q  +  a 


-  and  1  -  m  =  — — 

2q  +  a 


and 


1  +q 


c-J 

2q  +  2a 


c-j 


2q 


c-j 

q  +  a 


c-j 


2q  +  a 


c-j 


v2q  +  a 


c-j' 


(42) 

(43) 


q  -f  q 

=  (1  -  m)  - £zI.  (44) 

q 
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Using  (43)  and  (44),  (42)  becomes: 


u  =  (a  .  +  q)  cos 
c-J 


Let 


(1  -  n)(lp  -  0*)  +  -  Q  COS 


(1 


a  .  +  q 

M)-^J - 

q 


(<P  -  0*>  + 


P  -  (1  -  j u.)(<p  -  (//;|;), 

then  (45)  becomes 


u  -  (a  .  +  q)  cos()3  +  -  q  cos 

c  3 


Similarly  one  finds  that 


a  .  +  q 


(45) 

(46) 

(47) 


v  =  (a  _.  +  q)  sin(|3  +  i// J  -  q  sin 
c  3 


a  .  +  q 

- P  +  0  ^ 


(48) 


Equations  (47)  and  (4  8)  are  the  parametric  equations  for  an  epicycloid 
generated  by  a  circle  of  radius  q  rolling  on  a  circle  of  radius  a  .  (see 
appendix  7). 


Now  let 


$2  = 


then  from  (17): 


c-J 

q 


A  a  _  ' 

and  hence,  q  =  C  -- 

H 


(49) 


2  ,  2a 

q  + 


y  -  1 


A2 

=  q 


and 


A 

q  -  a 


c-J 


a 

=  c~.1 


c-J 


a  .  =  2q 

c-J 


(50) 


(51) 


From  (51)  it  is  seen  that  the  circle  of  radius  q  also  rolls  on  the  inside 
a  2  2  2 

of  the  circle  of  radius  q.  On  the  circle  u  +  v  =  a  _ the  speed  (q)  is  the 

C  J 

Chapman- Jouguet  sound  speed,  hence  this  circle  is  called  the  Msonic  circle, ” 

One  then  can  say  that  equations  (47)  and  (4  8)  are  epicycloids  generated  by  a 

a  i 

circle  rolling  between  the  sonic  circle  and  the  circle  of  radius  q  =  - 

A  ^ 

Note  also  that  from  (50),  when  q  =  q,  a  =  0  and  conversely.  This 

2  A2 

indicates  that  the  circle  q  =  q  is  the  locus  of  points  (x,  y)  where  there  is  no 
material  or  where  the  pressure  is  zero  (see  equation  (19),  appendix  1), 
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To  investigate  the  mapping  of  x,  y  points  to  u,  v  points  nondimensionalize 
(47)  and  (4  8)  be  letting 

(52) 


6  ,  A. 

u  =  a  .  77  and  v  =  a„  A 
c-j  2 


c-j  2 


and 


®2-?+1'  r3  =  ^-1 


(53) 


so  that  substituting  from  (40),  (41),  and  (44)  into  (47)  and  (48), 

'Rr 
£ 

r 

l3 

'Rr 
£ 

v~ 
l3 


6  =  R2  cos  Q3  +  <//*)  -  R3  cosf  —  13  +  0*), 

3 

X  =  R2  sin  ((3  +  0J  -  Rg  ®tn  (r~  ^  +  ^*)‘ 


(54) 


Using  (43)  and  (46), 

0  =  (1  -  /u)(ip  -  (//,;.) 


or 


(3 


(55) 


vl/2 


-  0*  =  r— a* 

From  (37)  one  concludes,  since  1  -  M  =  1  >  0,  that  rT  curves  are 

generated  by (3  >0,  and  F~  curves  are  generated  by  0  <  0. 

In  other  words,  r+  curves  are  generated  by  rolling  the  circle  of  radius 
q  counterclockwise  around  the  "sonic"  circle,  while  F  curves  are  generated 
by  clockwise  rotation.  0  is  the  angle  subtended  at  the  center  of  the  "sonic" 
circle  by  the  center  of  the  circle  of  radius  q  as  it  rolls  from  its  initial  posi¬ 
tion  to  its  current  position,  and  0^  is  the  angle  subtended  at  the  center  of  the 
"sonic"  circle  between  the  radius  vector  to  the  center  of  the  circle  of  radius 
q  at  its  initial  position  and  the  positive  u  axis  (see  appendix  7). 

From  (33)  and  (25), 


yP  -  0* 

and  from  (55), 


< 


w 


rH  ^  * 
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whence 

|0|  <h^J£>,  <56) 

and  for  y  -  3,  p  =  (2/ 4)^^  =  1/  so  that  ^  -  -  —  or 

2V2 

|  ft 1  <  «  37.26°. 

Suppose  one  considers  the  semi-infinite  slab  of  H.  E.  (shown  by  the 
dotted  line  in  Fig.  23a  before  detonation)  in  the  x,  y  plane.  In  appendix  4  it  is 

y 

shown  that  along  the  detonation  front  (x  =  0),  v  =  X  =  0  and  u  =  ac_j  “  ~  +  ^  D 
or  6  =2.  Hence,  the  entire  line  [1]  -[4]  in  the  x,  y  plane  is  mapped  to  the 
point  6  =  2,  X  =  0  in  the  6,X  plane  (see  Figs.  23a  and  b).  From  (54)  this 
corresponds  to  j3  =  0,  <//,.  =  0.  Along  y  =  0,  v  =  X  =  0  so  that  the  x  axis  is 

A 

mapped  to  the  u  axis  between  the  "sonic"  circle  and  the  circle  of  radius  q. 

Along  the  line  [4]  -[5]  ,  p  =  a  =  0,  and  as  stated  above  all  points  where  a  =  0 

are  mapped  to  the  circle  q  =  q  .  Since  there  must  be  contact  between  the 

lines  [4]  -[5]  and  [4]  -[1] ,  the  point  [4]  in  the  x,  y  plane  is  a  singular  point  and 

is  mapped  to  the  Tq  curve  that  runs  between  the  "sonic"  circle  and  the  circle 

of  radius  q  from  the  point  6  =  2,  X  =  0. 

To  see  that  [4]  -[5]  is  indeed  a  straight  line  in  the  x,  y  plane,  note  that 

(37)  implies  that  the  T+  and  T~  curves  are  generated  by  rolling  the  generating 

circle  counterclockwise  and  clockwise,  respectively,  so  that  the  T  curves  at 

2  A2 

the  point  [4,  5]  in  the  u,  v  plane  are  both  tangent  to  the  circle  q  =  q  there. 
Since  the  flow  direction  q  bisects  the  tangent  to  each  of  r  and  T  ,  q  is 
normal  to  this  circle  at  [4,  5] .  Referring  to  equation  (25),  appendix  2,  it  is 
seen  that  a  =  0  gives  the  slope  of  both  C  curves  since  v/u  =  tan  6  =  direction 
of  the  flow.  However,  the  flow  direction  bisects  the  angle  between  C  and  C 
curves  and,  in  this  case,  must  be  collinear  with  them,  so  that  the  image  in 
the  x,  y  plane  of  [4,  5]  is  the  flow'  direction  of  points  on  the  upper  boundary  and 
has  constant  slope. 

Considering  the  F^  curve  and  (56)  it  is  seen  that  the  value  of  ft  at  [4,  5] 
is  given  by 

0  _  7T (1  -  p) 

P  “  2p 


and  to  generate  all  other  T  curves  in  the  region  bounded  by  X  =  0,  Tn,  and 
o  ao  _l  r  _/i  ..\i  ./.+  u 


q^  =  q^  one  lets  ipt.  vary  from  0  to  — — J  =  v  *  with  the  lower  limit  of  ft 

for  each  value  being  determined  by  X  =  0.  Similarly,  the  T-  curves  are 


7T (1  -  p)]  _  ijjZ 


Fig.  23.  Mapping. 
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7r(l  -  U.)  //'- 

determined  by  letting  0;,.  vary  between  0  and  — - — —  =  with  /3  <  0  and  the 

limits  of  |3  determined  by  Tq  and  X  =  0  for  0  <  ip},.  <  jj- — —  and  by  T ^  and 

for  It  is  seen  that  points  in  the  6,  X  plane  are 

determined  by  assigning  a  pair  of  values  ipZ.,  ip I.  Note  also  that  along  X  =  0, 
the  intersection  of  a  r+  and  r~  curve  occurs  where  ipZ.  =  by  the  way  in 
which  these  curves  are  generated. 

Using  (43)  and  (46)  above: 

0 


if 


1  -  jU 


(57) 


so  that 


+ 


13 


+ 


+ 


.+ 


0.,,  and  ip 


+  if* 


M  r  *  r  1  -  M 

Hence,  by  using  (57)  and  (58)  in  (52)  and  (54),  one  can  write  equations 
(19)  of  appendix  5  as: 

ti+  2u 


(58) 


2v 

i 

c-J 


+  +x  /R  2  +  +\ 

-  r  R2  cos(j3  +i //.,,)  -  Rg  cosl^-  /3  +05;;) 
j  '  \  3  V 

7  =  R2  sin(/3+  +  ipt)  -  R3  sin  ^2.  jS+  +  ipt\ 


(59) 


C 


-1 


V 

2u 

3. 

C“J 

2v 


tan 


—  + 1 p... 
1-1 U 


x 


R2  cos (j3~  +  0l) 


R3  cos 


=  R2  sin(]3  +  4j Z)  -  R3  sint^  j3  + 
c-j  '  V  3  v 


(60) 


(61) 


(S' 


tan l — - — 
V1  -  p 


+  i/'t 

'  -i' 


X 


(S' 


(62) 


2  2  2 

In  appendix  2  (equation  (17))  it  was  stated  that  u  +  v  -  a  >0  was  the 
condition  for  a  hyperbolic  system  of  equations  and  this  condition  has  been 
assumed  when  (58)-(62)  are  given  as  the  solution.  However,  at  the  detonation 


front  itself  u  =  a 


c-J 


2  2  2 

0,  and  a  =  a  .,  whence  u  +  v  -  a 

c- J 


0  and  the 
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system  is  parabolic  with  the  front,  the  line  [1] -[4]  in  Fig.  23,  in  the  x,  y 
plane  being  both  a  C+  and  a  C~  characteristic.  This  means  that  the  solution 
in  the  region  of  interest,  i.  e.,  that  enclosed  by  the  lines  [2]  -[1] ,  [1]  -[4] , 

[4]  -[5]  in  the  x,  y  plane,  cannot  be  obtained  from  the  known  information  at  the 
front  since  no  characteristic  curves  intersect  the  front  and  move  into  the 
region  of  interest.  In  the  u,  v  plane  the  situation  is  different. 

Recall  that  the  detonation  front  in  the  x,  y  plane  is,  except  for  the  point 

[4] ,  mapped  to  the  single  point  u  =  a  .,  v  =  0  in  the  u,  v  plane  and  the  point 

[4]  is  mapped  to  the  line  FQ.  Therefore  along  rQ  the  values  of  x,  y  are 

known,  i.  e.,  x  =  0,  y  =  b.  Along  v  =  0  the  values  of  y  are  known  (y  =  0)  and 

along  q  =  q  ,  x  =  oo,  and  y  =  oo  while  at  u  =  a  .,  v  =  0,  x  =  0,  0  <  y  <  b.  If 

the  point  u  =  a  .,  v  =  0  is  avoided  then  the  problem  can  be  solved  by  working 
c- j 

in  both  planes. 

In  the  u,  v  plane  the  problem  is  completely  solved  for  u  and  v.  Near  the 
point  u  =  ac  .,  v  =  0,  the  F~  curves  connecting  FQ  and  v  =  0  are  short 
(relatively)  and  it  seems  clear  that  at  least  one  can  be  found  which  is  accurately 
approximated  by  a  straight  line.  The  values  of  j3  and  ip*  can  be  selected  for  a 
point  on  both  FQ  and  v  =  0  and  x  and  y  are  known  on  rQ  while  y  =  0  on  v  =  0. 

In  other  words  only  x  on  v  =  0  is  not  known.  However,  along  the  F  curve 
connecting  these  two  points,  equation  (62)  can  be  used  in  finite  difference  form 
to  find  that  x  value.  Graphically  this  corresponds  to  getting  the  slope  of  the 
C~  cxirve  corresponding  to  this  r  curve  and  extending  a  line  from  the  point 
[4]  in  the  x,  y  plane  with  this  slope  until  it  intersects  the  line  [1]  -[2]  (see 

+ 

Fig.  23c,  d).  With  this  value  of  x  on  [1]  -[2]  the  values  of  x  and  y  on  the  C 
curve  through  this  point  can  be  found  using  (62)  along  C  from  Tq  and  (60) 
along  C+  from  T-  curves.  In  Fig.  23d  the  angles  0.  and  02  are  measured 
between  the  normals  to  Tq  and  F^  and  the  negative  v  axis  respectively.  Their 
average  is  used  to  drop  the  line  from  [0,  b]  in  Fig.  23c  to  [x,  0]  thus 
determining  x^.  The  values  x2,  y2  are  determined  similarly  using  normals  to 
both  r“,  r~,  and  F^,  r|  as  indicated.  That  the  procedure  works  is  shown  in 
the  example  problem  considered. 
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APPENDIX  7 

DERIVATION  OF  THE  PARAMETRIC 
EQUATIONS  FOR  AN  EPICYCLOID 

An  epicycloid  is  defined  as  the  path  traced  out  by  a  point  on  the  cir¬ 
cumference  of  a  circle  rolling  on  another  circle.  To  determine  the  general 
equation  for  such  a  curve,  consider  the  situation  depicted  in  Fig.  24  which 
shows  a  part  of  the  path  of  the  point  Q  on  the  circle  of  radius  q  which  was 
originally  at  point  P  on  the  circle  of  radius  a  _ ..  The  two  circles  are  in 

c_j 

contact  at  point  R  and  since  there  is  no  slipping  allowed,  the  length  along  the 
inner  circle  from  P  to  R  (called  C)  is  the  same  as  the  length  along  the  other 
circle  from  Q  to  R. 


Fig.  24 


From  Fig.  24,  note  the  following  relations: 


0  = 


a 


C 

q 


whence 


a  ./3 
C“J_ 


(1) 


a  - 


q 


(2) 
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.  0  7r  a 

a  +  2y  -  tt  or  T  -  2  "  2  ’ 

(3) 

L  ORQ  =  a  +  y . 

(4) 

Combining  (3)  and  (4): 

L  ORQ  -  |  |  . 

(5) 

Using  the  law  of  sines  in  triangle  RSQ: 

R^  _  q 

sin  a  sin  y  * 

(6) 

Using  (3)  in  (6)  and  sin  29  =  2  sin  6  cos  9: 
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n  n  2  0 

Now,  2  sin  ^  cos  =  sin  d  and  2  sin  =  1  -  cos  0,  so  that  this  last 


expression  may  be  written  as 

a  .0 

r»  _ 


a  .0 


u'  =  a  .  cos  0  +  q  sin  0  sin  c-  J--  +  q  cos  0  -  q  cos  0  cos  -J 

C-J  n 


or 


u'  =  (a  .  +  q)  cos  0  -  q 
c-j 


a  .0 


q 

a  .0 


n  c-t  .  0  ■  c-j 

cos  0  cos - - - sin  0  sin  — — 


whence 


Similarly 


u'  =  (a  .  +  q)  cos  0  -  q  cos 
c-j  H 


q 

a  •  +  q 
c-j 


0 


a  •  +  q 

—  —  c-1 

v'  =  (a  .  +  q)  sin  0  -  q  sin - *- - 0  . 

c-3  n 


(12) 


(13) 


To  rotate  the  u'  ,  v'  axes  through  the  angle  -ip as  shown,  use  the 
relation 


u 


hjj 

+  iv  =  e  '  (u*  +  iv’  ),  where  i  =  v/^1. 


or 


u  =  u'  cos  -  v*  sin  ip},., 

v  =  u*  sin  ipt,.  +  v'  cos  ip i,.. 

Using  (12)  and  (13)  in  (15), 

ac_i  +  q 

u  r  (a  .  +  q)  cos  0  cos  ip  -  q  cos  - * -  cos  ipt,, 

c  J  q 

a  •  +  q 

-  (a  _  .  +  q)  sin  0  sin  < //.,,+  q  sin - f -  0  sin  ip 

c_j  ~ 


(14) 


(15) 


or . 


u  =  (ac_j  +  q)  cos  (0  +  "  q  cos 


a  •  +  q 

- 0  +  ip* 


(16) 


and  similarly 


/a  _ .  +  q  \ 

v  =  (ac  j  +  q)  sin  (0  +  -  q  sin^— -j - 0  +  i p^J 


(17) 
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It  is  clear  from  Fig.  24  that  /3  >  0  generates  curves  going  counterclock¬ 
wise  and  j3  <  0  generates  curves  going  clockwise,  and  that  to  sweep  out  all  the 
counterclockwise  curves  one  varies  i/j,,,  and  similarly  for  the  clockwise  curves. 
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APPENDIX  8 

THE  NEWTON -RAPHSON  METHOD" 

If  one  is  given  the  equation  f(x)  =  0,  and  lets  x1  +  h1  =  x  where  x1  is  the 


ith  approximation  of  x  and  h1  is  the  error  in  the_ith  approximation,  then 


h*  =  -  —  x  ■  I-  -  where  f '  (x)  =  ^ 


and 


f'  (xl) 


i+1  i  ,  ,  i  i  f(x1) 
x  =  x  +  h  =  x 


f  •  (x1) 


Let 


then 


f(/3)  =  sin 


f*  (j?)  =  cos 


0  +  ^-(0+  -  rjj  ) 


$  +  l)(ip+  -  0  ) 


R- 

wr 


sin 


R?  1  + 
rT^  +  #  > 


=  o. 


-  cos 


R 

LR3 


r-/3  +i(0+  -0~) 


so  that  to  solve  (3)  for  j3  one  solves  by  iteration: 


£i+1  ^  i3X 


sin 

/31  +  t|-(0+  -  0~) 

R3  . 

-  sin 

R2 

[rJ  P1  +  W+ 

cos 

[/31  +  i(0+  -  0  ) 

-  COS 

[r^1  +  ^+ 

si+1  -  e1 

until  - - t — <—  <  some  reasonable  limit. 


(1) 


(2) 


(3) 


(4) 


(5) 


See  ref.  8,  p. 


192. 
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APPENDIX  9 

BERNOULLI'  S  EQUATION  FOR  "WILKINS" 

EQUATION  OF  STATE 

The  first  law  of  thermodynamics  can  be  written  as: 

T  dr>  =  de  +  pd  (J-). 

For  isentropic  flow  dr)  -  0,  whence: 
dc  =  -pd(^). 

From  equation  (8),  appendix  1, 

de  +  u  du  +  v  dv  +  dp  +  p  d(^)  =  0, 
and  combining  (2)  and  (3): 

u  du  +  v  dv  =  -  dp. 

r 

Now  the  relation  (4)  holds  along  lines  where  dp  -  0,  i.  e.,  along  path 
or  flow  lines.  Hence,  (4)  may  be  written  as 

1  d  ,2  2v  _  1  dp  _  a2  dp 

2  dS  (U  +  V  >  ‘  "pHS  “35 


2  2  rp(S)  a2 

u  +  v  +  2  \  —  dp  =  constant, 

JP(SJ  p 

where  S,  S.,.  are  two  values  of  a  parameter  along  a  path  line, 

2 

a  -  sound  speed  defined  as  in  appendix  1. 

Equation  (5)  is  the  well-known  Bernoulli' s  equation. 

To  express  (5)  in  terms  of  the  Wilkins  equation  of  state,  recall  the 
definition  of  relative  volume, 

V  -  PQ/p, 

2 

and  the  expression  for  a  from  equation  (2),  section  II,  part  A: 


2  _  1  AQ 


0  Vv 


„2  -RV  <1+“,Cs 


-+  BRV  e"itv  + 
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Using  these  last  two  expressions,  it  is  clear  that  to  evaluate  (5)  the  following 
expression  must  be  integrated: 


5. 


p(S)  a2dp  „V 


p(SJ 


I 


V(SJ 


Ag+RBV  e-RV  + 

Vu+1 


dV. 


(6) 


The  integration  of  (6)  is  quite  straightforward  and  the  resulting  equation  is: 


„  *  2  ,  2  ,  2 
Const  =  u  +  v  + 


AQ 


P0  L  (Q  -  DV1 


1  RV  ^  + 

<JTT+B<V-i  «"RV  +  M 
K  uV 


(7) 
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APPENDIX  10 

SOLUTION  FOR  CHARACTERISTICS  IN  THE  HODOGRAPH  PLANE 
FOR  THE  "WILKINS"  EQUATION  OF  STATE 

For  the  plane  problem  using  the  Wilkins  equation  of  state  the  arguments 
in  appendix  6  are  unchanged  down  through  equation  (16).  As  indicated  in 
appendix  6  the  basic  equations  are  of  the  same  form  for  both  sets  of  curves 
(T)  and  equation  (16),  appendix  6,  can  be  written  as 


(1) 


where  ijj  is  the  angle  the  tangent  to  aF  curve  makes  with  the  positive  u  axis 
and  g  is  the  magnitude  of  the  velocity  component  normal  to  a  F  curve. 

If  equation  (7),  appendix  9,  is  evaluated  at  the  Chapman- Jouguet  point 
(detonation  front),  then 


q2  +  A 

P0 


AQ  Q  ■  +  B(V  +  l)e 
(Q  -  DV^1  R 


-RV  +  t1  + 


uv 


w 


=  D 


(2) 


where 


„  2,2 
D  =  a  .  +  — 
C~J  P0 


AO  /  is  -RV0  .  (1  +  ld)C 

Q  7T  +  B(Vc-i+E)e  + - - 


(Q  -  1)V 


W- 


c-3 


and 


2  2  2  2  2 
q  =  u  +  v  =  a  +  g  . 


Combining  (2),  (3),  and  (4): 


2  „  2  2 
g  =  D  -  a 


po  L  (Q  -  DVQ 

Differentiating  (5)  with  respect  to  ip: 


AQ -  4.  B  (v  +  4)  e-RV  +  (1  +  “>Cs 

-1  v  R' 

uV 


(3) 

(4) 

(5) 


dg  dV  _  (  da2  ,  2a2\  dV 

2gdV^7'V'^r  ~)^F  ' 

Using  (1): 

dg  dV  _ 
dV  "  a  * 


(6) 


Hence,  putting  (7)  in  (6): 


-97- 


whence 


dV 
d ip 


dip 


2ga 


^  2 
da 

2a2 

dV 

V 

1  da2 

2 

a 

2  dV 

V 

ga 

(8) 


(9) 


If  (9)  can  be  integrated,  then  ip  can  be  determined  as  a  function  of  V  and 
the  relations  (13)  and  (14)  of  appendix  6  can  be  used  to  determine  u  and  v,  and 
the  equations  for  C+  and  C  will  determine  x  and  y. 

-f  + 

From  the  considerations  in  appendix  6  it  is  seen  that  for  T  (and  ip  ), 
g  <  0;  and  for  T  (and  ip  ),  g  >  0.  If  the  symbol  g  is  considered  to  be  a 
positive  quantity,  then  letting 


F  (V ) 


,  ,  2  2 

1  da  a 

2  "dV"  “  V" 


(10) 


one  has 


+ 

* 


F(V)dV  +  <//* 


and 

ip~  =  ^  F(V)dV  +ipl. 


(ID 


(12) 


2 

Hence,  with  g  given  by  the  square  root  of  (5),  a  by  equation  (2),  section  II, 

part  A,  \jj  by  (11),  and  0  by  (12),  equations  (13)  and  (14)  from  appendix  6 

can  be  used  to  give: 

_4-  4-  .  4* 

I  :  u  =  a  cos  ifj  +  g  sin  ip  3 

+  +  (13) 

v  =  a  sin  ip  -  g  cos  ip  , 

r  :  u  =  a  cos  </y  -  g  sin  ip  , 

(14) 

v  =  a  sin  \p  +  g  cos  i p  , 

provided  the  integral  in  (11)  and  (12)  exists  and  is  finite. 

With  a  considerable  amount  of  work  it  can  be  shown  numerically  that 


-98- 


F(v)  < 


=  F'(v) 


where  p  <  1  and  M'  -  constant,  as  V 


For  reference: 


1  f  AQ(Q  +  I_)  +  BR2v2e-RV  +  (1  +‘j)(2Ju)Cs 

2  vQ  vw  1 

1  r  AO  2  -RV  (l+u)Cs]l/2 

J_  AQ  +  BRV^e  KV  + - — £ 

vF  VQ_1  v“ 


RV  ./ 2  _ 

M'  =  Me  C_J  ^~0/(^RB  Vc_.)  . 

Since  the  integral  of  F'  (v)  exists  for  p  <  1  for  an  arbitrary  upper  limit, 

the  integral  of  F(v)  would  also  exist  by  the  comparison  test  if  one  could  be 

sure  that  the  function  ip  was  single  valued  at  V  =  V  ..  That  this  is  indeed 

true  can  be  seen  by  considering  Fig.  5  of  UCRL-7797  which  indicates  that 

near  V  .  the  material  behaves  like  an  ideal  gas.  The  characteristics  for  the 
c-j 

ideal  gas  equation  of  state  were  shown  to  become  horizontal  at  V  =  V  .,  and 
the  same  result  should  follow  in  this  case. 

It  seems  clear,  therefore,  that  the  integral  in  (11)  and  (12)  has  a  finite 
value  and  can  be  evaluated  numerically.  This  point  is  discussed  further  in 
section  II,  part  B. 

At  the  intersection  of  a  T+  and  T~  curve,  the  relations  (13)  and  (14) 
yield  the  following  formulas  for  finding  ip  in  terms  of  ijj  : 


tan  ip 


2  2+  + 
(a  -  g  )  sin  ip  -  2ag  cos  ip 

2 - 2 - +  + 

(a  -  g  )  cos  ip  +  2ag  sin  ip 


Along  v  =  0,  (13)  and  (14)  yield: 


tan  ip 


To  find  ip~  and  ip  at  the  intersection  of  a  T  curve  and  v  =  0, 


a  sin  ip  +  g  cos  ip  =0 
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must  be  solved  for  ip~ .  This  is  done  by  the  use  of  the  Newton-Raphson  method 
described  in  appendix  8  to  find  V  from 


gtv1"1) 

a(V*_1)  tan  ip  (V*  1) . 


(18) 


where  V*  is  the  ith  estimate  of  V;  i p  and  ip+  are  then  determined  from  (16)  and 
(17). 

In  a  similar  manner  ip+  and  ip~  are  determined  at  the  intersection  of  T 
and  r-  by  finding  V  from: 


i-1 


fdV1-1)  [sin  i//+(V1~ 1 )  +  sin^V1"1)]  1 
\a(V*  ■*■)  [cos^^V*  ■*■)  -  cos  ip  (V*  *)] J 


and  ip 


+ 


and  ip  are  calculated  from  (11)  and  (12). 


(19) 


-100- 


REFERENCES 

1.  UCRL-67 97,  "Shock  Hydrodynamics,"  Wilkins,  M.  L.  ,  Feb.  1962. 

2.  Supersonic  Flow  in  Shock  Waves,  Courant  and  Friedricks,  Interscience 
Publishers,  Inc.,  New  York,  1948. 

3.  "An  Investigation,  by  the  Method  of  Characteristics,  of  the  Lateral 
Expansion  of  Gases  Behind  a  Detonating  Slab  of  Explosive,"  Hill,  R.  ,  and 
Pack,  D.  C.,  Proceedings  of  the  Royal  Society  of  London,  Vol.  194,  1947, 

pp.  524-541. 

4.  UCRL-7797,  "The  Equation  of  State  for  PBX  9404  and  LX-04-01,' 
Wilkins,  M.  L.  ,  Squire,  B.  ,  Halperin,  B.  ,  27  April  1 964. 

5  UCRL-7569-T,  "Write-up  of  Subroutine  ROMBRG,  a  General-Purpose 

Routine  for  Numerical  Integration,"  Fritsch,  Fred  N. ,  Oct.  1963. 

6.  Nonlinear  Theory  of  Continuous  Media,  Eringen,  A.  C. ,  McGraw-Hill 

Book  Company,  Inc.,  New  York,  1962. 

7.  Methods  of  Mathematical  Physics,  Vol.  II,  Courant,  R.  ,  Interscience 

Publishers,  Inc.,  New  York,  1962. 

8.  Numerical  Mathematical  Analysis,  Scarborough,  J.  B.  ,  Johns  Hopkins 


Press,  3rd  Ed.  ,  1955. 


-101- 


APPENDIX.  1 1 

TABLES  OF  OUTPUT  FOR  PRESSURE  CURVES 


Table  1.  Pressure  at  various  distances  from  the  centerline  of  the  H.  E.  as 
a  function  of  distance  from  the  front  in  units  of  the  initial  thickness  of  the  slab, 
for  the  ideal  gas  equation  of  state  in  plane  geometry. 


y  = 

o 

o 

!B 

y  = 

0.1 

y  = 

0.2 

y  = 

0.3 

y  = 

0.4 

X 

p 

X 

P 

X 

P 

X 

P 

X 

p 

0.0013 

0.3900 

0.0016 

0.3900 

0.0022 

0.3900 

0.0033 

0.3899 

0.0066 

0.3887 

.0189 

.3894 

.01  92 

.3893 

.0206 

.3886 

.0181 

.3877 

.0212 

.3777 

.0407 

.3873 

.0415 

.3866 

.0411 

.3846 

.0389 

,3794 

.0448 

.3422 

.0649 

.3831 

.0599 

.3830 

.0559 

.3800 

.0613 

.3655 

.0600 

.3140 

.0822 

.3790 

.0835 

.3765 

.0803 

.3705 

.0828 

.34  82 

.0744 

.2870 

.1127 

.3698 

.1055 

.3691 

.1113 

.3546 

.1035 

.3293 

.0916 

.2569 

.1299 

.3635 

.1138 

.36  59 

.1274 

.3452 

.1179 

.3157 

.1257 

.2084 

.1439 

.3579 

.1443 

.3527 

.1402 

.3373 

.1292 

.304  9 

.1361 

.1965 

.1666 

.3480 

.1662 

.3421 

.1510 

.3304 

.1606 

.2759 

.1447 

.1876 

.1853 

.3392 

.1839 

.3329 

.  1  904 

.304  5 

.1820 

.2573 

.1738 

.1624 

.2016 

.3311 

.1992 

.3247 

.2182 

.2861 

.1987 

.2437 

.2074 

.1411 

.2230 

.3199 

.2250 

.3103 

.2405 

.2717 

.2127 

.2331 

.2421 

.3097 

.2469 

.2979 

.2595 

.2598 

b 

.2459 

.1233 

.2596 

.3001 

.2662 

.2869 

.2763 

.24  96 

.2619 

.2006 

.2711 

.1144 

.2811 

.2882 

.2837 

.2768 

.2914 

.2407 

.2955 

.1822 

.2903 

.1086 

.3012 

.2771 

.2  999 

.2676 

.3052 

.2328 

.3060 

.1044 

.3201 

.2666 

.3224 

.2551 

.3181 

.2257 

.3219 

.1696 

.3427 

.2542 

.3431 

.2437 

.3415 

.2133 

.3441 

.1602 

.3600 

.2448 

.3626 

.2333 

.3627 

.2027 

.3634 

,1526 

.3604 

.0927 

.3811 

.2337 

.3811 

.2237 

.3821 

.1935 

.3808 

.1463 

.4017 

.2231 

.4046 

.2120 

.4002 

.1853 

.3966 

.1410 

.3967 

.0866 

.4219 

.2130 

.4215 

.2038 

.4172 

.1780 

.4248 

.1321 

.4250 

.0825 

.4419 

.2033 

.4433 

.1937 

.4413 

.16  82 

.4376 

.1284 

.4486 

.07  94 

.4617 

.1941 

.4592 

.1866 

.4638 

.1595 

.4614 

.1219 

.4692 

.0769 

.4814 

.1853 

.4799 

.1777 

.4852 

.1517 

.4831 

.1164 

.4876 

.0748 

.5012 

.1768 

.5002 

.1693 

.5058 

.1447 

.5033 

.1116 

.5044 

.0730 

.5209 

.1686 

.5201 

.1615 

.5190 

.1404 

.5223 

.1073 

.5199 

.0714 

.5408 

.1608 

.53  98 

.1541 

.  53  84 

.1343 

.5403 

.1035 

.5479 

.0686 

.5608 

.1533 

.5593 

.1472 

.56  34 

.1269 

.5576 

.1000 

.5608 

.0674 

.5810 

.1461 

.5835 

.1390 

.5816 

.1218 

.5822 

.0953 

.5849 

.0652 

.6015 

.1391 

.6027 

.1328 

.5996 

.1170 

.5979 

.0924 

.5962 

.0642 

.6222 

.1324 

.6220 

.1270 

.6172 

.1125 

.6207 

.0884 

.6177 

.0623 

.6390 

.1272 

.6412 

.1214 

.6403 

.1069 

.6426 

.0848 

.6379 

.0606 

.6604 

.1210 

.6605 

.1161 

.6575 

.1030 

.6639 

.0814 

.6571 

.0591 

.6821 

.1150 

.6799 

.1110 

.6802 

.0980 

.6846 

.0783 

.6845 

.0569 

.7043 

.1092 

.6993 

.1062 

.7027 

.0934 

.6982 

.0763 

.7019 

.0556 

.7223 

.1047 

.7239 

.1004 

.7194 

.0901 

.7183 

.0736 

.7187 

.0543 

.7407 

.1004 

.7437 

.0960 

.7417 

.0859 

.7380 

.0709 

.7432 

.0525 

.7595 

.0962 

.7637 

.0918 

.7585 

.0830 

.7639 

.0676 

.7590 

.0514 

.7834 

.0911 

.7  840 

.0878 

.7807 

.07  92 

.7830 

.0653 

.7822 

.04  97 

.8030 

.0872 

.7993 

.0848 

.8031 

.07*6 

.8020 

.0631 

.7973 

.04  87 

.8230 

.0834 

.8200 

.0811 

.8199 

.0731 

.8209 

.0610 

.8195 

.0472 

.8435 

.0797 

.8410 

.0775 

.8424 

.0698 

.83  97 

.0590 

.8413 

.0458 

.8591 

.0770 

.8622 

.0740 

.85  94 

.0675 

.8585 

.0571 

.8628 

.0444 

.8804 

.0735 

.8839 

.0707 

.8821 

.0645 

.8834 

.0546 

.8770 

.0434 

.9023 

.0702 

.9003 

.06  83 

.8994 

.0623 

.9020 

.0529 

.8980 

.0423 

.9246 

,0669 

.9226 

.0652 

.9225 

.0595 

.9207 

.0512 

.9189 

.0411 

.9418 

.0646 

.93  96 

.0630 

.94  00 

.0575 

.9395 

.04  95 

.9396 

.0399 

.9593 

.0623 

.9626 

.0601 

.9637 

.0550 

.9582 

.04  80 

.9602 

.0388 

0.9832 

0.0593 

0.9801 

0.0580 

0.9816 

0.0531 

0.9834 

0.0459 

0.9807 

0.0377 

y  -  distance  from,  centerline  in  units  of  initial  thickness, 
x  =  distance  from  front  in  units  of  initial  thickness, 
p  =  pressure  in  megabars. 

kfilank  spaces  were  left  where  the  calculated  data  did  not  fall  near  the  values  of  the  x  coordinate  being  used 
in  the  other  y  =  constant  lines. 
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Table  2.  Pressure  at  various  distances  from  the  centerline  of  the  H.  E. 
as  a  function  of  the  distance  from  the  front  in  units  of  the  initial  slab  thick¬ 
ness,  for  Wilkins  equation  of  state  in  plane  geometry. 


y  = 

o.oa 

© 

y  = 

0.2 

y 

=  0.3 

y  = 

0.4 

X 

P 

X 

P 

X 

P 

X 

P 

X 

P 

0.0047 

0.3904 

0.0064 

0.3903 

0.0088 

0.3901 

_b 

.0211 

.3893 

.0192 

.3894 

.0204 

.3886 

0.0222 

0.3857 

0.0265 

0.3656 

.0404 

.3863 

.0388 

.3860 

.0375 

.3  842 

.0393 

.3763 

.0436 

.3310 

.0601 

.3813 

.0593 

.3807 

.0618 

.3745 

.0637 

.3558 

.0592 

.2951 

.0788 

.3762 

.0779 

.3744 

.07  96 

.3649 

.0779 

.3411 

.0817 

.2462 

.1002 

.36  84 

.1013 

.3645 

.0981 

.3530 

.0977 

.3188 

.1021 

.2090 

.1216 

.3590 

.1208 

.3545 

.1212 

.3362 

.1162 

.2973 

.1250 

.1760 

.1417 

.3489 

.1400 

.3438 

.1396 

.3222 

.1417 

.26  84 

.1450 

.1537 

.1617 

.3378 

.1606 

.3311 

.1597 

.3061 

.1574 

.2515 

.1630 

.1376 

.1817 

.3257 

.1810 

.3178 

.1792 

.2904 

.1788 

.2303 

.1797 

.1254 

.2021 

.3128 

.2011 

.3042 

.2037 

.2708 

.2026 

.2091 

.1989 

.1138 

.2187 

.3019 

.2213 

.2  902 

.2229 

.2559 

.2214 

.1939 

.2198 

.1035 

.2399 

.2877 

.2416 

.2760 

.2418 

.2417 

.2396 

.1807 

.2390 

.0958 

.2618 

.2729 

.2622 

.2617 

.2605 

.22  82 

.2596 

.1677 

.2569 

.0897 

.2799 

.2607 

.2778 

.2509 

.2790 

.2154 

.2829 

.1542 

.2769 

.0839 

.2985 

.2481 

.2990 

.2366 

.2976 

.2033 

.2981 

.1462 

.2991 

.0784 

.3177 

.2354 

.3207 

.2224 

.3225 

.1880 

.3205 

.1356 

.3208 

.0738 

.3377 

.2225 

.3373 

.2119 

.3413 

.1772 

.3424 

.1263 

.3408 

.0701 

.3584 

.2095 

.3600 

.1980 

.3603 

.1670 

.3569 

.1207 

.3566  ' 

.0675 

.3801 

.1964 

.3775 

.1878 

.3795 

.1572 

.3785 

.1129 

.3824 

.0637 

.4027 

.1833 

.4015 

.1744 

.3989 

.1479 

.4000 

.1059 

.3993 

.0614 

.4205 

.1735 

.4201 

.1646 

.4188 

.1390 

.4216 

.0995 

.4159 

.0594 

.4389 

.1638 

.4392 

.1551 

.4390 

.1306 

.4433 

.0936 

.4402 

.0566 

.4580 

.1542 

.4590 

.1457 

.4596 

.1225 

.4579 

.0898 

.4561 

.0549 

.4780 

.1447 

.4794 

.1366 

.4808 

.1148 

.4  800 

.0846 

.4799 

.0525 

.4990 

.1354 

.5005 

.1277 

.5025 

.1074 

.5024 

.0796 

.5037 

.0503 

.5209 

.1262 

.5225 

.1191 

.5173 

.1027 

.5176 

.0765 

.5195 

.0489 

.5439 

.1172 

.5376 

.1136 

.5401 

.0959 

.5408 

.0720 

.5434 

.0469 

.5599 

.1114 

.5611 

.1055 

.5637 

.0  894 

.5565 

.0691 

.5594 

.0456 

.5765 

.1056 

.5773 

.1002 

.5798 

.0852 

.5806 

.06  50 

.5837 

.0437 

.6027 

.0972 

.6026 

.0926 

.5964 

.0812 

.5970 

.0624 

.6001 

.0424 

.6209 

.0918 

.6201 

.0877 

.6220 

.0754 

.6222 

.0586 

.6167 

.0412 

.6400 

.0864 

.6382 

.0829 

.6396 

.0716 

.6395 

.0562 

.6421 

.0395 

.6598 

.0812 

.6570 

.07  83 

.6577 

.06  80 

.6571 

.0538 

.6593 

.0383 

.6806 

.0762 

.6764 

.0739 

.6765 

.0645 

.6  843 

.0504 

.6768 

.0372 

.7023 

.0713 

.6965 

.0695 

.6957 

.0612 

.7030 

.04  82 

.7037 

.0355 

.7250 

.0666 

.7174 

.0654 

.7156 

.0579 

.7223 

.0461 

.7220 

.0344 

.7369 

.0643 

.73  92 

.0614 

.7361 

.0548 

.7420 

.0441 

.7408 

.0334 

.7615 

.0598 

.7619 

.0575 

.7573 

.0518 

.7623 

.0421 

.7600 

.0323 

.7  874 

.0555 

.7856 

.0538 

.7792 

.0489 

.7832 

.0402 

.7796 

.0313 

.8009 

.0535 

.7978 

.0520 

.8020 

.0461 

.8047 

.0383 

.7997 

.0302 

.8291 

.04  94 

.8231 

.0485 

.8257 

.0434 

.8157 

.0374 

.8204 

.02  92 

.8438 

.0475 

.8362 

.0468 

.8378 

.0421 

.8383 

.0356 

.8416 

.0282 

.8590 

.04  56 

.8634 

.0436 

.8629 

.0396 

.8617 

.0339 

.8634 

.0273 

.8746 

.0437 

.8776 

.0420 

.8759 

.03  84 

.8858 

.0323 

.8745 

.0268 

.9075 

.0401 

.906  9 

.0390 

.9026 

.0361 

.8982 

.0315 

.8973 

.0258 

.924  8 

.03  84 

.9222 

.0376 

.9164 

.0349 

.9237 

.0299 

.9208 

.0249 

.9427 

.0368 

.9379 

.0362 

.9449 

.0328 

.9368 

.02  91 

.94  50 

.0240 

.9612 

.0352 

.9541 

.0348 

.9597 

.0317 

.9638 

.0277 

.9574 

.0235 

0.9805 

0.0336 

0.9878 

0.0322 

0.9748 

0.0307 

0.9778 

0.0270 

0.9829 

0.0227 

ay  =  distance  from  centerline  in  units  of  initial  thickness, 
x  £  distance  from  front  in  units  of  initial  thickness, 
p  =  pressure  in  megabars. 

^Blank  spaces  were  left  where  the  calculated  data  did  not  fall  near  the  values  of  the  x  coordinate  being  used 
in  the  other  y  =  constant  lines. 


Table  3.  Pressure  at  various  distances  from  the  centerline  of  the  H.  E. 
as  a  function  of  the  distance  from  the  front  in  units  of  the  initial  thickness 
of  the  cylinder,  for  Wilkins  equation  of  state  in  cylindrical  geometry.  a 


y 

=  o.ob 

y  = 

0.1 

y 

-  0.2 

y  = 

0.3 

y  = 

0.4 

X 

F 

X 

P 

X 

P 

X 

P 

X 

P 

0.0502 

0.3840 

c 

.0519 

.3832 

0.0517 

0.3823 

0.0524 

0.3784 

0.0522 

0.3659 

0.0515 

0.3119 

.0539 

.3822 

.0540 

.3813 

.0560 

.3813 

.0563 

.3802 

.0554 

.3766 

.0565 

.3612 

.0584 

.3804 

.0587 

.3791 

.0584 

.3747 

.0609 

.3795 

.0611 

.3781 

.0615 

.3728 

.0610 

.3562 

.0637 

.37  82 

.0635 

.3771 

.0647 

.3707 

.0668 

.3770 

.0664 

.3758 

.0672 

.2724 

.0688 

.3747 

.0701 

.3756 

.0715 

.3734 

.0711 

.3663 

.0700 

.3458 

.0738 

.3740 

.0743 

.3640 

.0745 

.3404 

.0748 

.2549 

.0779 

.3722 

.0776 

.3704 

.0777 

.3616 

.0822 

.3706 

.0810 

.36  85 

.0847 

.3665 

.0844 

.3564 

.0837 

.3290 

.0868 

.3689 

.0877 

.3539 

.0884 

.3226 

.0923 

.3659 

.0932 

.3616 

.0919 

.3504 

.0903 

.2209 

.0986 

.3622 

.0981 

.3585 

.0989 

.3447 

.0978 

.3105 

.1060 

.3578 

.1034 

.3551 

.1030 

.3411 

.1026 

.3039 

0.1048 

0.1942 

.1144 

.3526 

.1091 

.3514 

.1121 

.3331 

.1122 

.2908 

.1173 

.3282 

.1170 

.2846 

.1241 

.3461 

.1224 

.3419 

.1289 

.3173 

.1269 

.2713 

.1356 

.3382 

.1309 

.3354 

.1356 

.3107 

.1368 

.2583 

.1406 

.3279 

.1424 

.2511 

.1493 

.3281 

.1518 

.3187 

.1505 

.2960 

.1532 

.2381 

.1658 

.3155 

.1588 

.2876 

.1596 

.2303 

.1798 

.2947 

.1801 

.2660 

.1737 

.2137 

.1861 

.2986 

.1932 

.2530 

.1898 

.1963 

.1978 

.27  85 

.2084 

.1778 

.2181 

.2767 

.2192 

.2595 

.2256 

.2216 

0.2184 

0.1685 

.24  54 

.2474 

.2461 

.2353 

.2454 

.2038 

.2915 

.2082 

.2988 

.1609 

.3236 

.1717 

0.3334 

0.1367 

0.3594 

0.1560 

.3845 

.1315 

0.4745 

0.0869 

aThis  table,  unlike  Tables  1  and  2,  includes  all  the  output  calculated  by  the  program, 
ky  =  distance  from  centerline  in  units  of  initial  thickness, 
x  =  distance  from  front  in  units  of  initial  thickness, 
p  =  pressure  in  megabars. 

cBlank  spaces  were  left  so  that  the  x  values  across  the  table  would  be  more  or  less  comparable. 
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This  report  was  prepared  as  an  account  of  Government 
sponsored  work.  Neither  the  United  States,  nor  the  Com¬ 
mission,  nor  any  person  acting  on  behalf  of  the  Commission: 

A.  Makes  any  warranty  or  representation,  expressed  or 
implied,  with  respect  to  the  accuracy,  completeness, 
or  usefulness  of  the  information  contained  in  this 
report,  or  that  the  use  of  any  information,  appa¬ 
ratus,  method,  or  process  disclosed  in  this  report 
may  not  infringe  privately  owned  rights;  or 

B.  Assumes  any  liabilities  with  respect  to  the  use  of, 
or  for  damages  resulting  from  the  use  of  any  infor¬ 
mation,  apparatus,  method,  or  process  disclosed  in 
this  report. 

As  used  in  the  above,  "person  acting  on  behalf  of  the 
Commission"  includes  any  employee  or  contractor  of  the  Com¬ 
mission,  or  employee  of  such  contractor,  to  the  extent  that 
such  employee  or  contractor  of  the  Commission,  or  employee 
of  such  contractor  prepares,  disseminates,  or  provides  access 
to,  any  information  pursuant  to  his  employment  or  contract 
with  the  Commission,  or  his  employment  with  such  contractor. 


